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show  in  particular  that  diffusion  driven  instability  will  occur  only  when  the  diffusion 
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We  then  show  some  numerical  simulations  of  pattern  formation  in  the  Gray-Scott  model 
using  MATLAB  programs  to  implement  the  Galerkin  Spectral  method. 
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Reaction-Diffusion  Equations  and  Morphogenesis  using  Galerkin  Methods 

USNA 

Trident  Project 
Benjamin  M.  Heineike 
Project  advisers: 

Professor  Reza  Malek-Madani 
Associate  Professor  Sonia  Garcia 

1.  Introduction:  Modeling  Morphogenesis 

In  the  context  of  embryology,  morphogenesis  is,  as  its  Greek  roots  suggests,  the  part  of 
the  process  of  development  where  shape  and  form  begin  to  emerge.  The  process  by  which  each 
nearly  identical  stem  cell  makes  the  initial  decision  to  take  the  unique  pathway  that  leads  it  to 
fulfill  its  role  as  a  specialized,  functional  cell  in  the  fully-grown  organism  is  still  largely  a 
mystery.  After  the  cells  have  taken  those  first  steps,  their  subsequent  development  is  better 
understood,  but  the  nature  of  the  initial  impulse  that  causes  one  cell  to  become  a  muscle  cell  and 
another  to  become  a  nerve  cell  is  just  beginning  to  be  understood. 

Two  years  before  Alan  Turing  died  in  1954,  he  published  a  paper  entitled  “The  Chemical 
Basis  of  Morphogenesis ”  in  the  Philisophical  Transactions  of  the  Royal  Society  of  London[  25] 
that  had  a  much  different  flavor  than  the  work  on  computing  machines  and  artificial  intelligence 
for  which  he  was  so  famous.  Today  his  paper  is  the  cornerstone  of  a  body  of  scientific  literature 
suggesting  that  the  basis  of  morphogenesis  lies  in  the  pattern  formation  capabilities  of  certain 
chemicals  thought  to  be  present  in  an  early  embryo,  which  Turing  dubbed  morphogens. 
According  to  the  theory,  by  the  time  the  genes  begin  to  create  proteins  in  a  growing  embryo, 
there  is  already  a  chemical  pattern  of  morphogens  present  in  the  background  of  the  tissue.  These 
morphogens  interact  physically  and  chemically  with  one  another  to  create  patterns.  These 
chemical  patterns  provide  initial  developmental  signals  to  the  genes  that  cannot  diffuse 
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themselves  between  neighboring  cells.  If  the  interactions  of  the  morphogens  were  effected 
primarily  by  chemical  reactions  and  physical  diffusion,  Turing  suggested  that  one  could  model 
their  pattern  fonnation  mathematically  using  physical  and  chemical  laws.  More  specifically,  one 
could  model  morphogenesis  with  a  particular  class  of  partial  differential  equations  called 
reaction-diffusion  equations.  In  1991,  Lewis  Wolpert  articulated  a  theory  of  positional 
information  in  his  book  The  Triumph  of  the  Embryo[27]  in  which  genes  were  influenced  by 
morphogen  concentrations  distributed  in  a  gradient.  In  one  example  the  location  of  each  digit  on 
a  chick’s  wing  was  thought  to  be  pinpointed  by  a  corresponding  concentration  of  the  morphogen, 
retinoic  acid,  which  decreased  as  it  diffused  away  from  a  source  along  the  antero-posterior  axis 
of  the  wing  (which  is  the  axis  determined  by  a  pencil  in  one’s  hand  if  it  were  held  with  a  closed 
fist). 

Only  recently  with  the  work  of  researchers  like  Ouyang,  Swinney,  and  Petrov  [17], [19] 
have  experimental  “Turing  Patterns”  been  exhibited  on  actual  biological  chemicals  in  a 
laboratory.  The  self-organizing  patterns  that  biologists  are  seeing  in  their  petri  dishes  and 
gelatinous  solutions  have  been  organizing  themselves  for  years  on  the  computer  screens  of 
mathematicians  [9], [14], [15], [18].  From  the  time  that  Turing  first  showed  that  nonhomogeneous 
solutions  existed  for  reaction-diffusion  systems  modeling  morphogenesis  in  a  simple  ring  of 
cells,  there  has  been  significant  interest  in  the  types  of  patterns  that  arise  in  the  solutions  of  these 
reaction-diffusion  equations. 

Reaction-diffusion  equations  are  partial  differential  equations  (PDEs)  which  model  the 
way  collections  of  particles  behave  taking  into  account  two  forces.  First,  they  model  the  way 
that  the  particles  interact  with  one  another  at  a  given  point  in  space,  described  in  chemical  terms 
as  their  reactions.  Second,  they  model  the  way  that  a  given  particle  distributes  itself  in  space, 
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diffusing  from  regions  of  higher  concentrations  towards  regions  of  lower  concentrations. 
Reaction-diffusion  equations  can  model  more  than  just  morphogenesis,  having  been  used  in  the 
past  to  model  population  densities,  host/parasite  models  [20],  electrical  reactions  that  occur 
between  nerve  cells  [8],  and  chemical  waves  such  as  those  found  in  the  Belousov-Zhabotinsky 
reaction  (ch.  7,  [15]).  Turing’s  remarkable  discovery  was  the  fact  that  the  presence  of  diffusion, 
which  acts  in  a  region  to  distribute  particles  more  homogeneously,  given  the  right  interactions 
between  particles  being  modeled,  could  help  initiate  very  interesting,  nonhomogeneous  patterns. 

Only  today  with  advanced  computer  technology  and  mathematical  methods  are  we 
beginning  appreciate  the  rich  pattern  fonnation  capabilities  of  these  models.  As  a  general  rule, 
the  variability  of  the  initial  and  boundary  conditions  as  well  as  the  complex  interaction  between 
concentrations  over  both  space  and  time  preclude  a  simple  algorithmic  solution  to  a  reaction- 
diffusion  system.  As  with  most  partial  differential  equations  used  in  modeling,  exact  solutions 
are  seldom  discovered  or  even  attempted.  Therefore,  to  investigate  our  reaction-diffusion 
system,  which  models  two  morphogens  in  two  dimensions,  we  will  turn  towards  approximate 
solutions.  Recent  advancements  in  scientific  computing  and  numerical  analysis  have  increased 
the  availability  of  high-powered  computing  resources  to  researchers,  rendering  approximate 
solutions  of  partial  differential  equations  an  efficient  and  reliable  tool. 

The  particular  approximation  method  that  we  will  use  is  called  the  Galerkin  Spectral 
method.  Many  approximation  methods  (the  Finite  Difference  method  for  instance)  break  the 
problem  into  smaller  workable  bits  by  dividing  the  spatial  domain  into  pieces  on  a  mesh  and  then 
putting  the  individual  puzzle  pieces  of  the  solution  back  together.  Instead  of  breaking  the 
solutions  down  into  puzzle  pieces,  spectral  methods  can  be  thought  of  as  breaking  the  solution 
into  layers.  We  assume  that  the  solution  to  the  system  can  be  written  as  a  linear  combination  of  a 
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product  of  unknown  functions  in  time  with  known  basis  functions  over  space.  The  Galerkin 
method  is  a  way  to  determine  these  unknown  functions  in  time  by  eliminating  spatial  dependence 
using  inner  products. 

Galerkin  methods  reduce  PDEs  to  a  collection  of  Ordinary  Differential  Equations  (ODEs) 
in  time.  For  an  exact  solution,  we  would  need  a  system  of  infinitely  many  ODEs,  but  for  our 
purposes  a  finite  approximation  will  be  sufficient  to  capture  the  overall  behavior  of  the  system. 
The  primary  advantage  of  the  Galerkin  method  is  that  it  is  relatively  easy  to  program  and  to 
increase  the  number  of  “layers”  in  the  approximate  solution.  As  will  become  clear  in  the 
subsequent  sections,  just  as  one  is  able  to  increase  the  accuracy  of  a  Fourier  series  representation 
of  a  solution  by  taking  more  terms  in  the  Fourier  basis,  we  are  able  to  add  more  accuracy  to  our 
approximate  solution  by  taking  more  terms  in  the  so-called  Galerkin  basis. 

The  focus  of  this  project  is  to  create  a  body  of  programs  in  MATLAB  that  can  solve  a 
system  of  nonlinear  reaction-diffusion  equations  using  the  Galerkin  method,  and  then  use  them 
to  analyze  the  pattern  formation  capabilities  of  a  particular  system  in  the  context  of 
morphogenesis.  In  particular,  we  are  looking  for  nonhomo geneous  steady-state  and  periodic 
solutions  of  the  Gray-Scott  model,  a  particular  two  species  reaction-diffusion  system  in  two 
dimensions  [18].  Along  the  way  we  use  the  method  to  test  certain  one-dimensional  nonlinear 
reaction-diffusion  equations  such  as  the  Burgers  equation  and  the  Allen-Cahn  equation  [4]. 

For  the  Gray-Scott  model,  we  will  also  analyze  the  ODEs  created  by  the  lower  degree 
Galerkin  method  solutions,  and  find  regions  in  parameter  space  near  which  particular 
equilibrium  points  of  the  solution  undergo  Hopf  bifurcation.  We  also  use  the  principles  of 
spectral  analysis  to  find  conditions  on  the  diffusion  parameters  that  exhibit  diffusion  driven 
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instability.  This  analysis  helps  us  to  find  parameter  values  for  which  the  solution  exhibits 
periodic  behavior  and  nonhomogeneous  steady-state  solutions. 

Intrinsic  to  this  study  is  an  awareness  of  the  complexity  of  these  systems.  In  a 
mathematical  context,  complex  does  not  simply  mean  complicated.  Rather,  complexity 
describes  the  phenomenon  whereby  complicated  structures  arise  from  the  simple  interactions  of 
interconnected  individual  particles.  The  structures  that  emerge  from  these  interactions  are  not 
virtually  random  like  chaotic  systems,  but  are  intricately  structured.  However,  the  eventual 
organization  of  a  system  will  not  immediately  follow  from  the  rules  of  interaction  for  the 
particles,  and  could  not  possibly  be  extrapolated  from  an  examination  of  a  particular  particle  by 
itself.  One  would  use  complexity  to  describe  the  way  that  the  interaction  of  each  individual  ant 
with  its  neighbors  and  the  environment  gives  rise  to  a  seemingly  conscious  colony,  or  how 
individual  water  and  air  molecules  acting  according  to  physical  and  chemical  principles  can  form 
highly  ordered  weather  structures  like  storm  clouds  and  tornadoes.  Our  morphogenesis  model  is 
yet  another  application  of  this  concept.  We  propose  that  some  of  the  complicated  patterns 
exhibited  in  living  organisms  with  such  a  high  level  of  reproducibility  are  controlled  by  the 
simple  interactions  of  morphogens  in  highly  interconnected  networks  of  individual  cells.  The 
purpose  of  studying  this  type  of  complexity  is  to  gain  clues  about  the  characteristics  of  an 
aggregate  structure  from  knowledge  of  the  way  individual  particles  interact,  without  resorting  to 
isolating  the  particle  from  its  neighbors.  We  would  like  to  know  not  only  how  an  individual 
particle  affects  the  entire  system,  but  also  what  characteristics  of  the  entire  system  cause  the 
particles  to  organize  themselves.  First  we  will  define  reaction-diffusion  equations  and  the 


model. 
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2.  Reaction-Diffusion  Equations 

A  system  of  reaction-diffusion  equations  is  a  system  of  equations  of  the  form 


—  =  DAu  +  f(u,Vu,x,0 
dt 


(2.1) 


over  a  region  DcR',  where  u(x,t)  is  a  vector  representing 
the  states  (in  our  model  morphogen  concentrations)  of  a 
group  of  substances  at  time  t  and  position  xeR".  A  is  a 
matrix  of  diffusion  coefficients,  which  in  a  two  species 


system  is  typically  of  the  form  D  = 


ju  0 
0  v 


,  and  Auis  the 


Laplacian  differential  operator  acting  on  u  with  respect  to 
xgD.  It  is  the  second  order  spatial  rate  of  change  of  u.  In  the  most  general  case,  the  inputs  in 
the  reaction  function  f  are  u,  Vu ,  the  gradient  of  u  with  respect  to  x,  x  and  t  (p.5,  [1]).  These 
partial  differential  equations  are  subject  to  boundary  conditions  over  Q  cz  D  and  initial 
conditions.  In  these  equations  the  term  containing  the  Laplacian  operator  is  the  diffusion  term. 
Without  the  function  f,  (2.1)  is  the  heat  equation,  one  of  the  first  equations  encountered  in  any 
partial  differential  equation  course.  The  heat  equation  models  the  diffusion  of  heat  from  regions 
of  higher  temperature,  or  heat  concentration,  to  regions  of  lower  temperature,  which  is  very 
similar  to  chemical  diffusion.  The  function  f  is  called  the  reaction  function  because  it  represents 
the  interactions  between  particles  that  act  to  increase  or  decrease  the  quantities  of  each  species, 
and  may  depend  on  the  concentration  of  particles  themselves  (u),  the  gradient  of  the 
concentrations  with  respect  to  space  ( Vu ),  and  the  location  of  the  reaction  in  space  and  time,  (x 
and  t).  The  use  of  chemical  terms  is  meant  merely  as  an  analogy,  as  reaction-diffusion  equations 
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have  found  broad  application  in  areas  other  than  chemistry,  such  as  neurological  signal 
transmission  [2],  Belousov-Zhabotinsky  chemical  waves  (pi 59,  [2]),  geochemical  systems 
(p229,  [2]  ),  combustion  theory  (pi  14  [2]),  and  other  complex  systems.  In  particular,  the 
Burgers  equation, 

ut  =  euxx  +uux , 

can  be  classified  as  a  reaction-diffusion  equation  and  is  one  of  the  most  important  equations  in 
the  study  of  fluid  dynamics.  We  will  use  it  later  to  test  our  numerical  methods. 

This  study  is  primarily  concerned  with  functions  f  that  depend  only  on  the  concentrations 
of  the  reactants.  This  idealization  is  a  good  approximation  for  many  chemical  reactions  held  at 
constant  temperature,  as  is  often  true  for  biological  reactions.  The  general  system  now  reduces 
to 

—  =  DAu  +  f(u).  (2.2) 

dt 


This  system  is  augmented  by  initial  conditions, 
u(x,0)  =  h(x) , 

and  boundary  conditions.  We  will  be  dealing  with  two  morphogens,  that  is,  the  vector  u  will  be 
given  by: 


u(x,f) 


u(x,t) 

v(x,0 


Murray  [11]  gives  a  good  overview  of  several  reaction-diffusion  equations  used  to  model 
morphogenesis.  The  three  primary  reactions  he  mentions  are  Schnakenberg’s  reaction,  Gierer 
and  Meinhardt’s  activator/inhibitor  model,  and  Thomas’  experimental  model.  Schnakenberg’s 
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model  is  mathematically  accessible,  but  has  not  found  much  biological  application.  It  is  given  in 
nondimensionalized  form  by 
u,  =  y(b  -v2u)  +  dAu, 

(2.6) 

vt  =  y(a  -  v  +  uv  )  +  Av. 

Gierer  and  Meinhardt  [10],  [5]  proposed  the  following  model  which  is  known  as  an 
activator/inhibitor  system. 

(  u2^ 

u, =y  a-bv-\ - +A u, 

{  v  J  (2.7) 

v,  =  y(u2  -  v)  +  dAv. 

The  nondimensionalized  parameters  are  the  same  as  above.  For  this  equation  we  will  call  u  the 
activator,  since  it  acts  to  increase  the  population  of  both  chemicals,  and  v  will  be  the  inhibitor, 
since  it  decreases  the  rate  of  change  over  time  for  each  morphogen.  For  patterns  to  occur,  Gierer 
and  Meinhardt  showed  that  d»  1,  in  other  words,  that  the  inhibitor  must  diffuse  significantly 
faster  than  the  activator.  As  an  illustration,  consider  a  predator/prey  system  [14].  Think  of  the 
prey  as  the  activators  and  the  predators  as  the  inhibitors.  The  predators,  cheetahs  for  instance, 
‘diffuse’  faster  than  the  prey,  antelope.  Where  the  antelope  gather  together  they  create  an 
environment  where  more  of  their  kind  can  thrive,  but  the  fast  moving  cheetahs  inhibit  their 
numbers  (through  digestion)  when  they  stray  from  the  herd.  Also  contact  between  antelope  and 
cheetahs  (again  through  digestion)  activates  the  production  of  cheetahs.  For  the  right 
parameters,  activator/inhibitor  reaction-diffusion  systems  form  dappled  patterns  where  activators 
clump  together  that  can  be  thought  of  as  analogous  to  herds  of  prey  species. 
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The  reaction-diffusion  system  that  we  will  focus  on  is  called  the  Gray-Scott  Model.  It 
has  been  studied  numerically  by  Pearson,  who  used  a  finite  difference  approach  [18].  The 
system  is  given  by 

u,  =  d„Au -uv2  +  F(l-u), 

(2.8) 

v,  =  dv Av  +  uv2  -  (F  +  k)v, 

where  k  is  the  dimensionless  rate  constant,  F  is  the  dimensionless  feed  rate,  and  d„  and  dv  are  the 
diffusion  coefficients.  For  our  simulations  they  will  typically  be  du= 2  x  10"5  and  dv=  10"5  which 
are  the  same  values  that  Pearson  used.  Our  goal  is  to  understand  the  effect  of  the  parameters  F, 
k,  du  and  dv  on  the  long  term  behavior  of  the  solutions  of  the  initial-boundary  value  problem 
associated  with  (2.8).  We  would  also  like  to  understand  the  effects  these  parameters  have  on  the 
accuracy  of  our  numerical  approximation  method.  More  concisely,  we  would  like  to  know  what 
causes  patterns  to  form  in  (2.8)  and  in  reaction-diffusion  equations  in  general. 
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3.  The  Galerkin  Spectral  Method 

Our  method  of  numerically  simulating  reaction-diffusion  equations  is  the  Galerkin  Spectral 
method.  Whereas  many  of  the  most  common  PDE  approximation  methods  used  today,  including 
the  finite  difference  method,  discretize  the  problem  in  space,  the  Galerkin  and  other  spectral 
methods  discretize  the  problem  over  a  spectrum  of  functions  that  are  continuous  over  the  whole 
space.  In  that  sense  spectral  methods  are  more  global  in  nature.  This  spectrum  of  functions, 
which  is  chosen  according  to  the  boundary  conditions,  forms  an  orthogonal  basis  for  the  function 
space  in  which  we  seek  the  solution  of  the  PDE.  To  implement  spectral  methods  we  construct  an 
approximate  solution  that  is  a  linear  combination  of  unknown  time  functions  coupled  with 
known  spatial  basis  functions  that  satisfy  the  boundary  conditions  associated  with  the  PDE.  The 
procedure  is  analogous  to  using  Fourier  analysis  to  express  any  bounded  continuous  function  as  a 
linear  combination  of  sine  and  cosine  functions  at  various  frequencies. 

The  concept  of  finding  a  set  of  basis  functions  with  which  we  can  represent  all  manner  of 
functions  that  satisfy  the  boundary  conditions  is  central  to  spectral  methods,  and  should  be 
explained  further.  In  finite  dimensional  vector  spaces,  a  basis  is  a  linearly  independent  set  of 
vectors  such  that  any  vector  in  the  vector  space  can  be  written  as  a  linear  combination  of  the 
basis  vectors.  To  illustrate,  let  Fbe  a  finite  dimensional  vector  space  over  the  field  A  and  B=(fi i, 
P 2....,  Pm)  be  a  basis  in  V.  Then  for  each  v  in  V  there  are  constants  a,-  e  A,  also  called  the 
coordinates  of  v  in  the  basis,  such  that 

v  =  a\P\  +  a1P1  + ...  +  amPm .  (3.1) 

Since  (3.1)  is  true  for  all  vectors  v  in  V then  we  say  that  B  spans  the  vector  space.  For  B  to  form 
a  basis  of  V,  it  must  do  more  than  just  span  the  entire  vector  space.  B  must  also  be  a  be  linearly 
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independent  set,  meaning  that  no  single  basis  vector  in  the  set  can  be  written  in  terms  of  the 
others.  In  that  sense,  a  basis  is  the  smallest  sized  set  that  can  span  a  vector  space.  We  often 
choose  basis  vectors  in  such  a  way  that  they  are  mutually  orthogonal,  that  is  if  (a,b)  represents 
the  inner  product  then 

(#,/?,)  =  0  when  i±j.  (3.2) 

Equation  (3.2)  will  help  us  compute  the  coordinates  a,  from  (3.1).  Taking  the  inner  product  of 
both  sides  of  (3.1)  with  a  fixed  basis  vector  /?,•  yields  (v,/f)  =  ai(j3i,fij )  or 


_  (y,Pi) 

(A.  A)' 


(3.3) 


An  important  property  of  a  set  of  basis  vectors  is  that  it  helps  define  what  is  called  a  weak 
characterization  of  the  zero  vector: 

v  =  0  if  and  only  if  (v,/7)  =  0  for  all  (3.4) 

We  will  see  later  that  a  similar  characterization  for  infinite  dimensional  function  spaces  will 
prove  essential  to  the  Galerkin  Spectral  method.  For  more  on  vector  spaces  and  basis  sets,  see 
Chapter  2  of  [7]. 

Finite  dimensional  vector  spaces  differ  from  their  infinite-dimensional  counterparts  in 
one  important  respect.  In  finite  dimensions,  as  long  as  the  number  of  linearly  independent 
vectors  in  a  set  is  equal  to  the  dimension  of  the  vector  space,  the  set  will  span  the  entire  space. 
The  fact  the  vectors  in  the  set  are  linearly  independent  and  span  the  vector  space  means  that  the 
set  is  a  basis  for  the  vector  space.  In  an  infinite  dimensional  vector  space,  however,  we  cannot 
be  sure  that  a  particular  set  is  a  basis  even  if  it  contains  infinitely  many  linearly  independent 
vectors  because  we  do  not  know  whether  or  not  the  vectors  span  the  space.  The  space  in  which 
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the  solutions  to  our  partial  differential  equations  lie,  the  set  of  bounded  functions  defined  over 
the  domain  Q  which  satisfy  the  boundary  conditions,  is  an  infinite  dimensional  function  space. 
We  must,  therefore,  be  thorough  when  detennining  whether  or  not  a  proposed  set  of  basis 
functions  spans  a  vector  space. 

Before  we  go  any  further,  let  us  define  orthogonality  for  our  infinite  dimensional  function 
space.  Let  J(x)  and  g(x)  be  real-valued  functions  in  the  space.  Their  inner  product  (fg)  is 
defined  as 

(f,g)  =  \f(x)g(x)dx.  (3.5) 

n 

When/ and  g  are  orthogonal  this  inner  product  is  zero. 

Let  C  =  {/j,^2,^3,...}  be  an  orthogonal  set  of  basis  functions  from  the  function  space  G. 

Since  the  functions  in  C  are  orthogonal,  they  automatically  satisfy  the  linear  independence 
condition.  In  addition  the  set  C  must  also  be  complete  in  G,  that  is,  it  must  also  be  possible  for 
any  function/in  G  to  be  expressed  as  a  linear  combination  of  this  infinite  set  of  basis  functions. 
More  concisely,  for  each /in  G, 

oo 

f(x)  =  ^a,//x)  (3.6) 

Z=1 

for  some  constants  {a\,  ai_,  03,...}. 

Equality  in  (3.6)  means: 

M 

lim  f(x)~Yjai(t>l(x)  =0, 

00  — - 

Z=1 

where 

\\gf  =(g,g)  =  \\gixf  dx. 

Q 
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Let  H  be  the  set  of  all  one-dimensional,  square  integrable,  real-valued  functions  defined 
on  [0,1]  that  vanish  at  the  endpoints: 

H  =  {f:R  ->  [0,1]  |  \\g\\  <  qo,  /(O)  =  /( 1)  =  0} . 

A  well-known  result  from  Fourier  analysis  states  that  the  set 

C  =  {sin(^x),sin(2^x),...,sin(n^x),...}  (3.7) 

forms  a  basis  for  H  ([6]  sec.  5.3.3).  Just  as  we  expressed  any  vector  in  a  vector  space  in  terms  of 
coordinates  with  respect  to  a  set  of  basis  vectors,  given  by  (3.3),  we  can  express  any  function  in 
H  in  terms  of  coordinates  with  respect  to  C.  Using  the  orthogonality  of  sin(/7rx)  over  the  domain 
[0,1],  which  means 

(sin (inx ) ,  sin (jnx ) )= 0  if  i^j,  (3.8) 

we  see  that 

=  WfMM  (3  9} 

Notice  that  if  we  remove  any  of  the  functions  in  C  with  an  even  coefficent  inside  the  sine 

function  to  fonn  a  new  set  C*,  we  will  still  have  infinitely  many  orthogonal  functions  that  satisfy 

the  boundary  conditions,  but  will  no  longer  have  a  basis  for  H  since  not  every  function  in  H  can 

be  written  as  a  linear  combination  of  the  functions  in  C*.  For  example,  consider  removing 

sin(27rx)  from  C  to  generate  C*.  If  C*  were  a  basis  for  H  then  we  could  write  sin(27ix)  in  terms 

of  its  coordinates.  However,  each  of  the  at  are  determined  by  (3.9): 

(sin(2;rx),sin(/;rx)) 

a,  = - =  0, 

(sin(z  n  x),  sin(z  n  x)) 
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where  i  can  be  every  natural  number  except  for  2.  Since  each  of  the  coordinates  is  0,  then 
sin(2roc)  would  have  to  be  the  0  function,  a  contradiction.  C*  cannot,  therefore,  be  a  basis  for  H 
even  though  it  consists  of  infinitely  many  linearly  independent  functions  in  H. 

The  Galerkin  Spectral  method  searches  for  solutions  of  the  PDE  system  in  terms  of  linear 
combinations  of  the  basis  functions  multiplied  by  unknown  functions  in  time.  The  fact  that  the 
domain  on  which  we  are  looking  for  the  solutions  is  bounded  assures  us  that  we  will  be  able  to 
form  a  basis  with  countably  many  basis  functions.  The  unknown  functions  in  time  are  then 
found  by  solving  the  now  countably  many  ODEs  that  are  created  by  substituting  the  template 
solution  into  the  PDE. 

For  the  reaction-diffusion  system  given  by  (2.2) 

—  =  DAu  +  /  (w) 

St  (3.10) 

we  define  the  differential  operator 

L(u)  =  ut  -  DAu  -  f(u).  (3.11) 

Note  that  (3.10)  is  equivalent  to 

L[»]=0.  (3.12) 

We  search  for  solutions  to  (3.10)  of  the  form: 

oo 

u(x,t)  =  YJai(t)</>i(x),  (3-13) 

1=1 

where  the  functions  satisfy  the  appropriate  boundary  conditions.  To  find  the  unknown  a,(t), 
we  now  substitute  (3.13)  into  (3.12)  and  use  the  characterization  of  zero  similar  to  (3.4): 
g(x)  s  0  if  and  only  if  (g(x),^(x))  =  0  for  all  i.  (3.14) 


18 


This  characterization  is  computationally  useful  because  in  order  to  show  that  a  function  is  zero 
for  each  value  in  its  continuous  domain,  one  must  only  show  that  the  countably  many  (but 
infinite)  conditions  in  (3.14)  are  satisfied,  namely  that  the  inner  product  of  each  basis  function 
with  g  is  zero.  Setting  Z(w )  =  0  in  this  sense,  we  are  left  with  the  countable  conditions: 

(^,Z[m])  =  0, 

(^2,Z[m])  =  0,  (3.15) 

The  inner  product  operation  in  (3.15)  removes  spatial  dependence  and  these  conditions  leave  us 
with  a  system  of  ODEs  in  terms  of  the  unknown  coefficients  at{i).  In  practice,  instead  of  using 
infinitely  many  basis  functions  in  (3.13),  we  truncate  the  solution  template  at  some  finite  number 
N.  Thus  in  (3.15)  we  would  have  N  ODEs  in  N  unkowns.  To  find  the  initial  conditions  needed 
to  solve  (3.15),  we  use  the  initial  conditions  in  the  original  PDE, 
u(x,0 )=h(x)  . 

Since  (3.13)  must  satisfy  these  initial  conditions  as  well,  we  impose  the  following  conditions  on 
the 

oo 

2>,(0#,(jc)  =  w(x,  0)  =  w(x,0)  =  h{x). 

i= 1 

Thus  the  initial  condition  for  each  function  <2/(0,  61,(0),  is  given  by  its  coordinate  for  h(x)  in  terms 
of  the  basis  functions.  By  (3.9)  we  have  an  explicit  formula; 

iU  (4(jc),*,(jc))' 

Now  that  we  have  N  ODEs  in  N  unknowns  with  initial  conditions,  we  can  solve  for  the 


a/(0  and  then  use  them  to  construct  an  approximate  solution  for  the  PDE.  We  can  solve  these 
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ODEs  numerically  using  any  number  of  ODE  solving  packages  found  in  mathematical  software 
such  as  MATLAB  and  Mathematica. 

We  will  now  illustrate  the  steps  associated  with  the  implementation  of  the  Galerkin 
method  on  the  following  basic  reaction-diffusion  equation  in  one  space  dimension: 

ut  =duxx+f(u),  (3.16) 


where  d  is  a  constant  real  number.  The  spatial  domain  is  D=[0,1].  The  equation  in  (3.16)  is 
subjected  to  Dirichlet  boundary  conditions, 

u(0,t)  =  u(l,t)  =  0,  (3.17) 


and  initial  conditions, 
u(x, 0)  =  g(x)  . 


(3.18) 


Step  1.  Choose  the  basis  functions  and  solution  template:  The  primary  precondition  for  this 
choice  is  that  the  basis  functions  satisfy  the  boundary  conditions.  If  C  =  is  our 

basis  then  the  solution  template  is  of  the  form: 


N 

u(t,x)  =  YJan(t)<t>n(x). 

n= 1 


(3.19) 


For  the  boundary  conditions  in  this  example  we  will  choose 


(/>„(x)  =  sm(n7rx). 


(3.20) 
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where  n  is  an  integer  since  (3.20)  is  zero  when  x=\  and  x=0.  Other  considerations  for  this  choice 
have  to  do  with  the  type  of  solution  that  one  expects,  the  information  about  the  solution  that  we 
are  trying  to  extract,  and  of  course  the  level  of  complexity  that  our  numerical  simulator  is 
prepared  to  handle.  Certain  basis  functions  may  be  more  natural  for  the  solution;  a  continuous 
solution  would  probably  be  better  represented  in  a  continuous  basis  set  of  functions.  The 
individual  modes  of  the  basis  functions  may  have  physical  interpretations  for  the  model,  and  we 
can  choose  a  basis  set  that  takes  advantage  of  these  interpretations  to  reveal  certain  properties 
about  the  model.  Later,  when  we  look  at  the  Gray-Scott  model,  the  0  mode  solution  will 
represent  the  system  without  diffusion,  and  we  can  ‘add’  diffusion  to  our  approximations  by 
taking  the  solution  with  more  modes.  Table  (3.1)  gives  examples  of  Fourier  basis  functions  one 
can  use  for  the  boundary  conditions  listed  in  the  left  column. 


Table  3.1 


Type  of  Boundary  Conditions 

Fourier  basis  to  use 

Periodic 

innx 

e 

Dirichlet:  u(0,t)=0,  u(l,t)=0 

sin(nnx) 

Neumann:  ux(0,t)=0,  ux(  1  ,t)=0 

cos(htdc) 

Mixed:  ux(0,t)=0,  u(l,t)=0 

cos 

f  (1  +  2  n)7Dc^ 

V  2  ) 
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There  are  a  variety  of  basis  functions  that  will  work  with  a  given  set  of  boundary  conditions.  If 
we  did  not  want  to  use  Fourier  basis  functions  for  one  reason  or  another,  we  could  look  for 
Chebyshev  and  Legendre  polynomial  bases  as  well  as  wavelets. 

Step  2.  Substitute  template  into  differential  operator  and  obtain  a  set  of  ODEs  in  time: 

Our  differential  operator  will  be 

L(u)  =  u,  -  duxx  -/(«)  .  (3.21) 


Substituting  u  from  (3.19)  into  (3.21)  we  have: 

N  f  N  A  f  N  \ 

L(u)  =  ^a'n(t)sm(n7ix)  +  d  ^n27r2an(t)sm(n7ix)  -f  ^a„(t)sin(n;zx)  .  (3.22) 

v»=i  J  v»=i  J 


n= 1 


2  2 

Notice  that  the  second  derivative  has  been  replaced  in  (3.22)  with  the  constant  -n  n  because 
(sin(«7rx))iv  =  -n2 n2  sin(/?7r x) .  In  order  to  obtain  the  a„(t),  we  set  (3.22)  equal  to  zero  using 
the  weak  formulation  described  in  (3.14)  and  (3.15): 


^ a '  (t)(sin{mnx),sin(n7Dc))+  7i2 d^ci" {t)n2 (sin(m^x),sin(/7^x))-  sin (rrmx),f  \  (t)sin(n®c) 

n=l  V.  \n=l  J) 


n= 1 


(  N 


=  0 


m 


=  {l,2,...,M- 


(3.23) 


Notice  that  the  inner  product  operation  and  the  summation  operations  commute  because 
integration  is  a  linear  operation.  The  functions  (ftjx)  =  sin(/?/zx)are  orthogonal  on  the  interval 
[0,1]  and  when  m=n,  (sin(/H/rx),  sin(/?7rx))=l/2  so  (3.23)  now  takes  the  fonn: 
a'm  (t)  -n2m2d 


<t>(x)m,f  Yjan^ni.x) 

V H=1  J  J 


(3.24) 
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for 


Step  3.  Find  the  Initial  conditions  for  the  ODEs:  In  order  to  find  a  unique  solution  to  a  system 
of  ordinary  differential  equations,  we  need  initial  conditions  for  each  function  in  the  system.  For 
the  initial  conditions,  we  use  (3.9)  to  get 


«.(0)  = 


(sin(m^x),g(x)) 
(sin  (mm),  sin(  myzw )) 


=  2(sin  (mnx),g(x)) 


(3.25) 


for 


Step  4.  Solve  the  initial  value  problem  for  time:  Now  we  have  N  ODEs  in  N  unknowns  with 
initial  conditions.  If  the  function/is  nonlinear,  then  the  system  of  ODEs  will  be  nonlinear,  and 
we  will  be  forced  in  most  cases  to  look  for  an  approximate  solution  using  a  numerical  ODE 
solver  such  as  MATLAB’s  ODE45  or  Mathematica’s  NDSolve.  The  main  computational  costs 
will  arise  from  the  last  tenn  in  (3.24)  which  is  the  nonlinear  term.  We  use  forward  time  stepping 
methods  primarily,  such  as  the  Runge-Kutta  methods  built  into  the  MATLAB  and  Mathematica 
ODE  solvers.  Many  reaction-diffusion  equations  have  a  property  called  stiffness  that  causes 
forward  time  stepping  methods  to  become  unstable  unless  the  time  step  is  made  extremely  small. 
MATLAB  and  Mathematica  use  highly  adaptive  routines  that  can  usually  avoid  this  problem; 
however,  the  computational  cost  of  solving  stiff  differential  equations  could  become  prohibitive 
even  for  adaptive  solvers. 

Step  5.  Reconstruct  the  solution:  Once  we  have  computed  all  the  coordinate  functions  an(t), 
we  can  reconstruct  the  approximate  solution  using  (3.19).  Since  we  solved  the  ODEs 
computationally,  our  a„(t)  will  not  be  continuous  functions  but  rather  discretized  tables  of  an(t) 
values  given  at  time  intervals  specified  by  our  numerical  ODE  solving  scheme.  Errors  in  our 
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final  approximate  solution  will  arise  at  several  points  during  the  process.  First  and  foremost, 
they  will  be  due  to  choosing  N  too  small.  The  magnitude  of  N  that  one  must  use  to  bound  the 
error  at  a  given  value  for  any  particular  reaction-diffusion  equation  is  a  difficult  problem  in 
analysis  that  will  not  be  addressed  here.  Errors  will  also  arise  from  numerically  approximating 
solutions  to  the  ODEs.  When  comparing  the  approximate  solution  to  actual  reaction-diffusion 
systems  found  in  nature,  we  must  also  take  into  account  errors  that  arise  from  approximations 
and  simplifications  that  were  made  when  creating  the  mathematical  model. 

The  tacit  assumption  made  throughout  this  study  is  that  our  initial-boundary  value 
problem  has  a  unique  solution.  This  assumption  can  be  made  rigorous  in  the  setting  of  some  of 
our  problems  (for  example,  the  nonlinear  Allen-Cahn  equation)  but  remains  a  formidable  open 
problem  in  mathematical  analysis  for  others  including  the  Burgers  equation  and  the  Gray-Scott 
model.  Indeed  the  complex  structure  of  the  periodic  and  steady-state  numerical  solutions  we 
obtain  is  an  indication  of  how  difficult  it  will  be  to  develop  an  existence  and  uniqueness  theory 
for  nonlinear  PDEs. 

Example  1:  One-dimensional  heat  equation 

The  simplest  reaction-diffusion  equation  is  one  that  does  not  contain  any  reaction  at  all.  This 
equation  is  commonly  referred  to  as  the  heat  equation  since  it  has  been  used  extensively  to  model 
the  way  heat  spreads  in  various  media  over  time.  It  is  given  in  one  dimension  by  (3.16)  where 
f(u)= 0: 


ut-kuxx ,  xe(0  ,L) 


(3.26) 


which  can  be  expressed  in  operator  notation  as: 
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L(u)  =  u,  -  ku ^  =  0  (3.27) 

with  initial  conditions 

u(x,  0)=g(x).  (3.28) 

In  (3.26)  k  represents  the  thermal  diffusivity  constant.  We  consider  here  Dirichlet  boundary 
conditions 


u(0,t )  =  u(L,t)  =  0 


(3.29) 


and  follow  the  step  by  step  process  outlined  above. 

Step  1-  Choose  the  basis  functions  and  solution  template:  Since  we  have  Dirichlet  boundary 

conditions,  we  choose  </>n  =  sin(-^)  .  Our  solution  template  is: 

N  n 

u(t,x)  =  '^jan(t)  sin( — x).  (3.30) 

n= 1  L 

Step  2-  Substitute  template  into  differential  operator  and  obtain  a  set  of  ODEs  in  time: 
Substituting  (3.30)  into  (3.27)  we  have 


Zf / \  *  / 

an(t)  sm(— 

n= 1  L 


x)  +  k 


I 

n= 1 


\ 

,  .  .  ,nn  . 

a»  (0sm(—  x) 


0. 


Using  the  weak  characterization  of  0  given  in  (3.14)  and  (3.15),  this  expression  gives  us  the  N 
conditions: 


N 


N 


2X(0  sin(— x),sin(— x)  +k^a„(0| 

n= 1  V  E  L  J  „=1 


^  H7T^ 


V  L  J 


V 


.  ,mn  .  .  ,nn  . 

sin(— —  x),  sin(-^—  x) 


=  0.  (3.31) 
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Using  the  fact  that: 


f 

f inn  ^ 

( nn  ^ 

)  1 

\L/ 2 

m  =  n 

sin 

- x 

,sin 

—  x 

V 

V  L  j 

V  L  J 

J  1 

[o 

m  ^  n 

(3.32) 


(3.31)  is  reduced  to 


<(0  =  ~k 


inn 

V  L  j 


a„M) 


since  the  constant  ^  factors  out  of  either  side. 


(3.33) 


Step  3-  Find  the  Initial  conditions  of  the  ODEs:  Using  (3.9)  and  the  orthogonality  of  our  basis 
set,  we  have: 


«„(0)  = 


sin 

V  L _ j_ 

f  .  ,mn  .  .  ,mn  f 

sin( - x),sin( - x) 

V  L  L  j 


=  Y^n(~-x)g(x)dx .  (3.34) 


The  <7,„(0)  are  nothing  more  than  the  Fourier  sine  coefficients  of  g(x) . 


Step  4-  Solve  the  initial  value  problem  for  time: 


Each  of  the  ODEs  in  (3.33)  has  the  solution 


am(t)  =  am(  0)e 


-k{'ny/Lf  i 


which  can  be  verified  by  direct  substitution. 

Step  5-  Reconstruct  the  solution:  Using  our  solution  template,  (3.30),  we  have 
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—  /  \  V"1  /rv\  I  .  H7T  . 

u(t,x)  =  2_jan(0)e  A;  sin(— x), 

«=i  -6 

where  a„(0)  is  the  nth  Fourier  sine  coefficient  of  the  initial  condition.  This  is  indeed  the  .Y-th 
partial  sum  of  the  true  solution  of  the  heat  equation  with  Dirichlet  boundary  conditions  which 
can  be  derived  by  more  traditional  methods  such  as  separation  of  variables. 

Example  2:  The  Burgers  equation  -  A  one-dimensional  nonlinear  reaction-diffusion  equation 

We  will  now  demonstrate  how  the  implementation  changes  when  a  nonlinear  reaction  is  added  to 
the  problem.  A  common  nonlinear  reaction-diffusion  equation  used  in  the  modeling  of 
turbulence  and  airflow  is  the  Burgers  equation, 

L  (u  )  =  u  t  -  £ u  xx  +  uu  x  =  0  .  (3.35) 

In  the  Burgers  equation  the  quantity  u  is  related  to  the  density  of  the  fluid  being  modeled  and  the 
term  s  uxx  represents  viscous  dissipation  in  the  model.  It  is  known  (see  Malek-Madani  [13], 

418-427)  that  when  s  =  0  that  the  typical  initial  disturbances  in  u  fonn  discontinuities  known  as 
shock  waves  and  yet  when  s>0  solutions  that  start  out  smooth  remain  smooth  for  all  time.  The 
computational  challenge  in  this  problem  is  to  gain  insight  into  the  Galerkin  method’s  response 
for  small  but  positive  viscosities.  We  will  analyze  (3.35)  over  the  interval  [0,1]  subject  to 
Dirichlet  boundary  conditions  on  the  same  domain  as  the  heat  equation  above  (3.29)  with  L=l. 
The  initial  condition  we  consider  is: 


|l-cos(8;zx)  re  [1/4, 1/2] 
u(x,0)  =  g(x)  =  ] 

0  otherwise. 


(3.36) 
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Just  as  before,  the  boundary  conditions  lead  us  to  choose  the  basis  (j)j  =  sin(//zx) .  Proceeding 

through  the  Galerkin  method  using  steps  1-5  and  the  relationship  (3.32),  we  arrive  at  the 
equations: 

a'  (t\  —  n2m2 f  (  N  N  3 

~jnz~  = - - - am (0 -  </>,„ OX S' a„ (Mn 0)Z aj Wj O) 

£  ^  \  n— 1  7=1  y 

We  can  simplify  the  nonlinear  term  for  ease  of  programming: 

</>„,  ox  S  a«  o)S  a'j  (Wj  o)  = Z  Z  a«  o  K-  (ofcx,  ox  k  0¥y  o)) = 

y  «=1  7=1  J  n= 1  7=1 

where  a  is  a  1  x.A/  row  vector  whose  entries  are  the  ah  i=l,..,N,  and  P  is  the  matrix  defined  by 

P,n  (X ./)  =  (f  „  OX#  0)#  O) ) 


It  is  easy  to  show  that 


7  +  m  =  U 

\j  ~  m\  =  /, 

(7  +  tn  &  i )  (  7  _  ^  /). 


Therefore,  the  Galerkin  method  reduces  to  N  differential  equations  in  N  unknowns  of  the  form 


a'Jt)  =  -7r2m2£am(t)-2(aPma') . 


The  initial  conditions  for  the  above  system  are  given  by  (3.9)  as 

am  (0)  =  (sin(m^x),  gO))  =  2  r  sin(jn7a^x  _  cos(8 nx)dx. 
(sin(/H/zx),sin(/tt/zx)]  ]-l4 


One  of  the  main  questions  we  have  regarding  the  Galerkin  Spectral  method  has  to  do  with  the 
magnitude  of  N  that  we  must  pick  to  accurately  describe  the  system.  This  choice  will  depend  on 
many  parameters,  the  most  important  being  the  amount  of  information  we  need  to  extract  from 
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our  approximation.  Figure  3.2  shows  Galerkin  approximations  to  the  Burgers  equation  for  8,  32 
and  64  basis  functions  at  time  1.  We  see  that  as  N  is  increased,  the  approximate  solution 
becomes  smoother  and  displays  less  oscillations  due  to  the  basis  functions. 


Burgers  Equation,  ut=eps*uxx  -uux 
eps=.001 ,  at  time  1 ,  u(x,0)=1-cos(8*Pi*x)  on  [1/4 ,1/2],  0  elsewhere 


Figure  3.2  Galerkin  approximation  to  the  Burgers  equation  for  various  N,  £=.001,  t=l. 

Figures  3. 3-3. 5  show  Galerkin  Spectral  method  approximations  to  the  Burgers  equation  for 
various  N  and  e  values.  In  these  diagrams,  each  curve  represents  our  approximate  solution  at 
time  intervals  of  0.05  units  with  a  final  time  of  1  unit.  With  N=32  and  s=0.01,  the  solution 
seems  to  qualitatively  follow  the  behavior  that  we  would  expect  from  the  Burgers  equation  (Fig 
3.3).  There  are  only  slight  oscillations  near  the  base  of  the  cosine  function  in  the  initial  condition 
for  the  first  time  value  due  to  the  inability  of  our  smooth  basis  functions  to  accurately 
approximate  the  sharp  corners  of  the  non-smooth  initial  condition.  In  Figure  3.4  where  we 
decrease  s  to  0.001,  and  keep  N  at  32,  the  Galerkin  method  shows  much  more  oscillations  and 
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inaccuracies  than  it  did  with  a  =  0.01  in  Fig  3.3. 
However,  when  we  increase  N  to  96,  we  see  that 
the  solution  begins  to  behave  more  smoothly 
despite  some  lingering  instabilities.  When  a  is 
increased,  the  solution  at  time  one  has  a  sharper 
comer  and  the  peak  of  the  solution  seems  to 
decrease  at  a  slower  rate.  This  is  because  a  is 
inversely  proportional  to  the  Reynold’s  number. 


Caleikin  Approximation  to  Burgers'  equalion 
sine  basis,  N=32,  eps=  001 

■jpi.l)  u(1.l>  0. 11(1,0)  1  cosjB*pi*i)  on  |1/4.1JV|  II  ebwtwui 


Each  Contour  ime 
rufiriisuritu  it  turm  iitnp 
of  OS  from  0  to  I 


Figure  3.3  The  Burgers  equation  n=32,  £=.01 

An  a  value  of  zero  would  model  a  fluid  with  no 


viscosity  (inviscid  flow),  and  a  shock  wave  would 
form.  With  a  positive  a  value,  the  shock  dissipates. 
Notice  also  that  the  oscillations,  which  are  caused 
by  the  inability  of  our  basis  functions  to  adapt  to 
non-smooth  shock  wave  behavior  begins  to  form 
for  smaller  a  values  decreases  with  time  as  a  result 


Figure  3.4  Burgers’  equation  n=32,  £=.001 


of  the  diffusion  term.  As  s  is  decreased,  the  effect 
of  the  diffusion  term  is  decreased,  so  we  expect 
oscillations  to  occur  earlier  and  last  longer.  This 
evidence  clearly  indicates  that  our  choice  of  N 
depends  on  the  parameter  s  in  the  equation. 


Figure  3.5  Burgers’  equation  n=96,  £=.001 
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Example  3:  The  Allen-Cahn  Equation 

This  next  example  illustrates  the  use  of  nonhomogeneous  boundary  conditions  as  well  as 
the  difficulty  of  knowing  when  one  has  reached  a  steady-state  solution  for  some  reaction- 
diffusion  systems.  The  Allen-Cahn  Equation  is  given  by 

L  (u  )  =  u  t  -  su  xx  -  u  (1  -  u  2  )  =  0  .  (3.37) 

Referring  to  (3.10)  we  havc/(w)=w(l-w").  Just  as  Trefethen  does  in  [23],  we  will  consider  the 
nonhomogeneous  Dirichlet  boundary  conditions: 

u  (  -  1 ,  t )  =  -1,  u  (1 ,  t )  =  1  ,  (3.38) 

and  the  initial  condition 

u(x, 0)  =  g(x )  =  0.53  x  +  0.47  sin(  - x) .  (3.39) 

Note  that  the  initial  condition  satisfies  the  boundary  conditions.  Additionally,  g(x)  is  asymmetric 
about  the  origin. 

Step  1-  Choose  the  basis  functions  and  solution  template: 

With  nonhomogeneous  Dirichlet  boundary  conditions,  we  choose  the  solution  template  in  the 
following  form: 

N 

u(t,x)  =  x  +  y^fln(Q  sin  {nn  x).  (3.40) 

n= 1 

Adding  the  term  x  to  the  solution  template  ensures  that  it  satisfies  the  boundary  conditions 
(3.38).  This  choice  is  by  no  means  unique,  so  we  have  added x,  the  first  function  that  came  to 
mind  that  satisfied  the  conditions.  The  basis  functions  will  still  be  r/;i;  (x)  =  sin  {nn  x)  but  notice 
that  the  solution  template  is  not  made  up  entirely  of  linear  combinations  of  the  basis  functions. 
Step  2-  Substitute  the  template  into  the  differential  operator  and  obtain  a  set  of  ODEs  in  time: 
Perfonning  the  procedures  used  in  (3.22)  and  (3.23)  we  are  left  with  ODEs  of  the  form: 
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(  N 


Yja'n  (0(^»  0)>  <t>n  (*))  +  n ' an  (0»  2  ($L  OX  <!>n  O))  +  ^  m  (*X  /  *  +  X  ttn  (0  tn  0) 

«=1  n=l  fy  \  n=l  y  J 


=  0 


m=l,...JV. 


(3-41) 


Since  the  basis  functions  are  orthogonal,  (3.41)  reduces  to: 


a'n  (t)  =  -7z1dm1am  (t)  + 


<t>,„ OX/  *  +  Ea»(o^„o) 


(3.42) 


v  v  »=i  yy 

The  last  term  in  (3.42)  is  the  most  complicated  term  and  will  be  the  most  costly  to  evaluate 
numerically  since/is  a  nonlinear  function. 

Step  3-  Find  the  Initial  conditions  of  the  ODEs:  Paralleling  the  procedure  we  used  to  get  (3.9), 
we  must  now  take  into  consideration  the  function  x  that  we  added  to  satisfy  the  boundary 
conditions.  We  start  with  the  initial  condition, 


N 

u{x,  0)  =  g(x )  =  x  +  ^  an  (0)  sin  (nnx), 

n= 1 


and  take  the  inner  product  with  each  basis  function: 

N 

^  a n  (0)(sin(/;z/zx),  sin(«/zx))  =  (sin {mm),  g(x)  -  x). 

n= 1 


Since  the  basis  functions  are  orthonormal,  meaning  that  they  are  orthogonal  and  (</>n ,  (j)n )  =  1  for 
all  n,  our  initial  conditions  are: 

a„,  (0)  =  (sin(m;zx),[g(x)  -  x]) .  (3.43) 

Steps  4  and  5  -  Solve  the  initial  value  problem  and  reconstruct  the  solution 

Now  as  above,  we  can  solve  the  initial  value  problem  (3.42)-(3.43)  in  MATLAB  (see  appendix). 


Figures  3.6  and  3.7  show  the  approximate  solution  using  the  Galerkin  method  with  8  basis 
functions.  These  figures  exhibit  an  important  phenomenon  known  as  bistability  found  in  some 
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reaction-diffusion  systems.  Notice  in  Figure  3.7  that  before  time  45,  the  solution  seems  to  have 
settled  into  a  nonhomogeneous,  wave-like  steady-state  solution,  but  suddenly  snaps  into  its  final 
nonhomogeneous  steady-state  solution  where  u  transitions  once  from  -1  to  1.  This  property  of 
reaction-diffusion  equations  makes  it  difficult  to  say  using  only  computational  evidence  when  a 
solution  has  reached  its  steady  state  and  when  it  is  merely  lingering  in  some  intermediate  state. 

Allen-Cahn  Equation 
sine  basis,  n=8 
1.5 

1 

0.5 

§  0 

-0.5 

-1 

-1.5, 

Figure  3.6  Allen-Cahn  Equation  Figure  3.7  Allen-Cahn  Equation 

Example  4:  The  heat  equation  in  two  dimensions 

A  few  changes  are  required  to  reformulate  our  problem  for  two  dimensions.  We  consider  the 
heat  equation  in  two  dimensions  with  the  heat  diffusivity  constant,  k  in  (3.27)  equal  to  one: 

L(u(x,y,t))  =  ut  -  Au  =  0  (x, y)  e  Q  =  [-1,1]  x  [-1,1], t  e  R+  (3.44) 

with  initial  conditions: 

u(x,yS))=g(x,y)  (3.45) 

and  periodic  boundary  conditions.  The  solution  template  takes  the  form: 

oo  oo 

U(x,y,t)=  Y,  Yamn^yt>,nn  k)  ' 

m—  co  n=- co 


(3.46) 
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As  in  one  dimension,  the  basis  functions  <f>„m(x,y )  are  determined  by  the  boundary  conditions.  To 
satisfy  periodic  boundary  conditions  we  choose: 

A,M.y)  =  e'1"""’'1.  (3.47) 

Notice  that 

(ftmnix-.y)  =  (/sin(/«/zx)  +  cos(  in  jix  ))(/s  i  n  ( /?  ny )  +  cos(nTry)) 
and 

A <t>mn(x,y)  =  -m27r2e,{m,a+nv) -n27r2e,(mmc+nv)  =  -{ni2n2  +n2n2)<f>mn{x,y)  .  (3.48) 

2  2  2 

In  other  words,  (j)mn  is  an  eigenfunction  of  the  Laplacian  operator  with  eigenvalues  -(m  +  n)n  . 
Substituting  (3.46)  into  (3.44)  and  using  (3.48),  we  have: 

oo  oo  oo  oo 

L{U)=  Y  Y^'^mn  (*W)  +  Y  Y Qmn  ^  Mnn  (*>  f)  •  (3-49) 

m=— co  n=— co  m=— oo  «=— oo 

We  now  take  the  inner  product  of  (3.49)  over  Q  with  (f>mn .  In  two  dimensions,  where  complex 
valued  functions  are  allowed,  the  inner  product  is  defined  by 
(f(x,y),g(x,y))  =  \[f(x  ,y)g(x,y)dA, 

Q 

where  g(x,y)  represents  the  complex  conjugate  of  the  function  g(x,y).  Approximating  the 
infinite  sums  with  finite  sums  between  -N  and  N,  we  will  be  left  with  (2 N+l)  ordinary 
differential  equations  in  t  of  the  form: 
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iii  a'mn  0VmH  0,  V),  (j)pq  (X,  V)l  =  I  I  ±a  mn  (t)n 2  (- m 2  -  n 2  )<f>mn  (x,  y),  </>pq  (x,  y ) 

\m=—N  n=—N  J  \m=—N  n=—N 


or  equivalently: 


N  N 


N  N 


X  Z«l(0(^™,(^v),^(x,v))=  X  X«,™(0^2(-^2  -n2)(<l>mn(x,y),<t>pq(x,  v)). 


m=—N  n=—N 


m=—N  n=—N 


(3.50) 


We  see  that  the  inner  product, 


t i>mn(x’y)’<t>Pq(x’y ))’ 


2 

is  a  known  quantity  for  any  given  p,q,m,  and  n  and  thus  we  have  a  system  of  (2N+  If  ODEs  in 

2 

(2 N+l)  unknowns  (the  amn(t )).  For  our  region  Q,  the  following  relationship  holds: 


k,,(x,yU„(x,y))=  }  }  e^+’^ei(^dxdy= 


0, 


(m  =  p)and(n=q), 
(m  ^  p)or(n  ^  q). 


(3.51) 


Using  (3.51)  we  can  rewrite  (3.50)  as 


a'mAt)  =  amn(t)X2(-™2  ~  n~)  , 


(3.52) 


m,n  e  {—N,...,N}  .  To  find  the  initial  conditions  with  orthogonal  basis  functions  we  use  the 
same  process  as  in  the  one-dimensional  case.  Evaluating  (3.46)  at  t= 0,  and  equating  the  initial 
condition  of  the  solution  template  to  that  of  the  true  solution  given  by  (3.45),  we  have: 
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X  2a™  =  i7(x’T’°)  =  W(X’T>°)  =  gU  v) .  (3.53) 

m=— oo  «=— oo 

We  take  the  inner  product  of  (3.53)  with  each  basis  function,  and  for  each  ^  every  one  of  the 
apq{ 0)  disappears  except  <7 ,/()),  resulting  in 

an  (°)  =  ^  (g(x,  y),  </>v  (x,  y)) .  (3.54) 

With  a  system  of  ODEs  and  initial  conditions,  we  again  have  an  initial  value  problem  for  which 
we  can  obtain  a  numerical  solution.  This  system  is  much  larger  than  the  one  our  scheme 
produced  for  the  one-dimensional  equations,  and  increasing  the  number  of  basis  functions  will 
take  significantly  more  computer  time  and  memory.  Even  for  the  reactionless  heat  equation  it 
takes  a  few  minutes  to  produce  an  8  basis  approximation  for  a  solution  at  time  100.  On  the  other 
hand,  the  system  of  ODEs  (3.52)  with  initial  conditions  (3.54)  could  be  solved  analytically 
because  (3.52)  is  an  uncoupled  system.  Once  we  add  our  nonlinear  reaction  terms,  it  will 
become  impossible  to  uncouple  the  ODEs  produced  by  the  Galerkin  method,  and  we  will 
normally  be  required  to  use  numerical  methods  to  solve  the  system. 

Example  5:  The  Gray-Scott  Model\ 

This  system  mentioned  in  (2.8)  will  be  the  focus  of  our  search  for  patterns.  It  is  a  variant  of  the 
autocatalytic  Selkov  model  of  glycolysis,  due  to  Gray  and  Scott([5],[21]). 


ut  =  du  Au  -  mv9  +  F(l  -  u ) 
v,  =  dvAv  +  uv2  -  (F  +  k)v 


(3.55) 
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In  the  model,  u  and  v  represent  the  concentrations  of  chemicals  U  and  V,  in  standard  chemistry 
notation  u=[U]  and  v=[V],  The  rate  of  change  of  the  concentrations  with  respect  to  time  is 
determined  both  by  diffusion,  modeled  by  the  Laplacian  of  the  concentrations  with  respect  to 
space,  and  the  chemical  reaction  rates  which  depend  on  the  concentrations  of  the  other  chemicals 
at  any  point.  The  diffusion  coefficients,  given  by  du  and  dv,  are  multiplied  by  the  diffusion  terms 
for  each  concentration.  The  reactions  involved  in  the  Gray-Scott  model  can  be  simplified  to  the 
following  reactions: 

U+2V^3V  (3.56) 

V^P.  (3.57) 

Both  reactions  proceed  in  only  one  direction.  Reaction  (3.56)  proceeds  at  a  rate  proportional  to 

2  2 

ravr  =uv  and  acts  to  decrease  the  concentration  of  chemical  U  and  increase  the  concentration 
of  chemical  V.  Reaction  (3.57)  converts  V  to  the  inert  product  P  at  a  rate  of  k\V\.  F  is  a  non- 
dimensionalized  feed  rate.  We  have  four  parameters  for  this  system,  F  and  k,  the  reaction 
parameters,  and  du  and  dv,  the  diffusion  coefficients. 

In  1993,  John  Pearson  from  Los  Alamos  National  Laboratories  published  results  from 
finite  difference  simulations  of  the  Gray-Scott  model  [18].  He  used  a  mesh  of  256  by  256  grid 
points  and  simulated  solutions  for  200,000  time  steps.  With  the  powerful  supercomputers  at  Los 
Alamos  National  Laboratories,  he  was  able  to  produce  some  very  interesting  patterns,  from 
oscillating  labyrinthine  stripes,  to  stable  hexagonal  patterns  of  points,  to  dividing  chemical  rings 
resembling  dividing  cells.  We  have  found  some  evidence  of  these  remarkable  patterns  using  the 
Galerkin  method  on  much  smaller  computers  for  much  smaller  amounts  of  time. 
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As  in  Pearson’s  simulations  we  start  with  diffusion  coefficients  of  du= 2  x  10'5  and 
<iv=10'5,  but  unlike  his  1993  paper,  we  show  patterns  that  occur  with  different  diffusion 
coefficients  as  well.  Our  simulations  take  place  on  the  region  Q=[-l,l]x[-l,l]  with  periodic 
boundary  conditions.  We  choose  the  initial  condition  as  a  small  perturbation  from  the 
homogeneous  equilibrium  point  (w0,vo)=(l,0).  Unlike  Pearson’s  step  function  perturbations,  our 
perturbations  will  be  smooth.  Explicitly,  the  perturbations  will  be  Gaussian  spots  (a  depression 
from  1  for  u  and  a  raised  impression  with  a  base  at  0  for  v).  The  following  is  an  example  of  an 
initial  condition  with  a  maximum  perturbation  magnitude  of  1/16  and  centered  at  ( x,y)=(h,k )  for 
both  chemical  u  and  chemical  v: 

16  (3.58) 

v(x,v,0)  =  0  +  -e-20<(^)2+(^)2). 

16 

A  homogeneous  equilibrium  point  is  a  point  at  which  the  concentrations  are  equal  throughout  the 
domain  and  are  not  changing  with  time.  Later  we  explain  how  we  know  (1,0)  is  a  homogeneous 
equilibrium  point  and  we  show  that  it  is  linearly  stable  for  all  reaction  and  diffusion  parameter 
values.  For  now  it  suffices  to  know  that  at  a  stable  equilibrium  point  the  concentrations  of  all  the 
chemicals  in  the  system  will  remain  constant  over  the  entire  domain,  and  will  return  to  that 
constant  state  when  slightly  perturbed  until  some  outside  force  acts  to  push  the  concentrations  far 
enough  away  from  their  equilibrium  values. 

The  purpose  of  using  initial  conditions  that  are  small  perturbations  from  homogeneous 
equilibrium  points  is  to  mimic  the  embryo’s  transition  from  a  stable,  patternless,  homogenous 
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state  to  a  state  where  stable,  nonhomogeneous  patterns  exist.  One  of  the  great  breakthroughs  of 
Turing’s  work  in  morphogenesis  was  that  he  showed  examples  of  homogeneous  steady-state 
solutions  that  appear  stable  when  only  the  chemical  reactions  are  modeled,  but  become  unstable 
when  diffusion  is  considered.  The  system  could  then  move  to  another  steady-state  solution, 
which  may  be  nonhomogeneous  [25].  Thus  a  chemical  system  initially  in  a  state  of 
homogeneous  equilibrium  could  take  on  a  variety  of  patterns,  dependant  on  the  reaction  and 
diffusion  parameters,  as  a  small  perturbation  of  the  initial  homogeneous  state  propagates 
throughout  the  domain.  Murray  explains  this  phenomenon  of  diffusion  driven  (or  Turing) 
instability  found  in  reaction-diffusion  equations  more  thoroughly  in  chapter  14  of  his  text 
Mathematical  Biology  [15]. 


Our  solution  template  is,  as  in  the  heat  equation: 

oo  oo 

u(x,y,t)=  Yj  Yja™(ty>mn(x’y)  (x,y)eQ  =  [-l,l]x[-l,l],teR+ 


oo  oo 

v(x,y,t )  =  E  'ZK.mjx.y),  (3.59) 

m=- oo  /7 =—oo 


where  (j)nm  (x,y)  is  given  by  (3.47).  We  will  demonstrate  the  Galerkin  procedure  for  the  first 

equaition  in  (3.55).  The  procedure  is  identical  for  the  second  equation.  Substituting  (3.59)  into 
the  linear  operator  defined  by  the  first  equation  in  (3.55)  we  have 


oo  oo 


oo  oo 


L(u )  =  EE  (x,y)  +  du  E  E“  mn(t)7l2(m2  +  n2)<f>m„(x,y) 


m=- oo  n=-  oo 


m=— oo  n =—oo 


00  00 


oo  oo 


+  EE»  mn^mn^yi  X  X Kn  W™  (*>  f) 

m=-co  «=— co  \ in  1= — co  nl=—<x>  J 


f  co  oo  \ 

-f  i-  E  Ea  mni^mni^y) 

V  m=— co  n=- oo  / 


=  0. 
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Using  the  weak  characterization  of  0  by  setting  the  inner  product  with  each  basis  function  (j)pq 
equal  to  zero,  we  are  left  with  equations  of  the  form: 


I  Z«-  (0(^„„,(uv),^w(x,y))  = 

m=-  oo  n=-  co 


00  00 


<  z 1”  m„  (t)7T 2  (m 2  +  n 2  )(</>mn  (x,  y),  (f)pq  (x,  y))  - 

m=-  oo  n=-  co 

oo  oo  /  oo  oo 

Z  I*  ml/?l  W  mini  (u  v)  ,<t>pq(x,y) 


m=— co  n=—co 


\m  1=— oo  nl=— oo 


+ 


00  00 

(F^pq(x,y))-f7'Z  ™(0(L(^v),^?(x,y)) 

m=—co  n=— oo 


(3.60) 


The  most  complicated  part  of  this  expression  contains  the  nonlinear  tenns  that  come  from 
substituting  the  template  solution  into  wv  .  Simplifying  this  expression,  as  we  did  with  the  heat 
equation,  we  have: 


(uv2,</>pq(x,y)))  = 


oo  oo 


oo  oo 


II«  run  (UT)I  X  ,(f>pq(x,V) 


m=-co  72 =—oo 


\  m =—oo  72 =-oo 


00  00  00  00  00  00 


Z  Z  Z  Z  Z  Z  ^  772  1  72 1  ^  722  2  72  2  ^  772  3  72  3  (*  7721/21  mini  (*,y)<i>  mini  (x,y),<l>pq(x,y)). 


772 1 = — OO  72 1= — 00  7722=— 00  722=— 00  7723=— 00  723=— 00 


Note  the  inner  product  inside  the  nonlinear  tenn  simplifies  to: 

(*>  y)<t>m2n2  (*>  TM„3„3  (*>  0,  k))  = 


I4’ 

0, 


(ml + m2 + m3 +p-0)  and  (n\ +n2+n3+q  =  0) , 
(ml + m2 + m3  +  ^  0)  or  (nl + m2 + m3 + q  ^  0) . 


(3.61) 
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If  we  truncate  the  infinite  sums  keeping  only  N  terms,  we  will  be  left  with  (2/V+  l  )2  equations  in 


2 

(2A/+1)  unknowns.  Denoting  the  sum 


N  N  N  N  N  N  N 

X  X  X  X  X  X  b>'  X-we  can  write 


m\=-N  n\=-N  m2=-N  n2=-N  m3=-N  n3=~N  M(N) 


(3.60)  in  terms  of  ODEs  of  the  form: 


00  00 

El  a'mn  (0(^™  (*>  y),  <!>„  (x,  y))  = 

m=—  oo  n=-co 
00  00 

X  Xa™WK,^2(-w2 

m=-co  n=- co 

wl«l  m2n2  (*,y)</>  m3n3  (.x,y),4>pq(,x,y))  + 

M(N) 

(F,(t>pq(x,y)) 


(3.62) 


Using  similar  reasoning,  we  can  write  the  second  equation  in  (3.55)  as: 


X  X  *»»  (0(^™  (u  V),  </>pq  (x,  y))  = 

m=—N  n=—N 

X  Xam„(0K^2(-^2  -n2)-F-%,„(x,v),f)w(x,v))+  (3.63) 

m=—N  n=-N 

WmW  (U  y)<f>m2„2  (x,  V)#  mM  (x,  V),  (j)  pq  (x,  y)). 

M(N) 

In  this  form  (3.62)  and  (3.63),  given  appropriate  initial  data,  can  be  solved  numerically  on  a 
computer  package  such  as  MATLAB.  This  is  made  easier  by  using  (3.61)  which  relies  on  the 
orthogonality  of  the  basis  functions  to  calculate  the  large  inner  product.  We  can  also  calculate 
the  inner  products  ahead  of  time  and  store  them  to  be  referenced  for  particular  index  values. 
Since  the  time  derivative  is  alone  on  the  left  side,  we  can  use  MATLAB ’s  built  in  ODE  solvers 
to  calculate  the  values  at  a  given  time.  See  the  appendix  for  the  programs  used  in  calculating  the 


results. 
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4.  The  Gray-Scott  Model 

The  Galerkin  Spectral  method  and  our  MATLAB  programs  provide  us  with  a  way  to  produce 
numerical  solutions  of  the  Gray-Scott  model  for  given  parameter  values.  We  will  now  analyze 
the  model  to  determine  where  in  parameter  space  we  might  begin  looking  for  patterns.  The  first 
step  is  to  examine  its  homogeneous  equilibria.  Recall  that  the  Gray-Scott  model  is  given  by: 

u.  =  d„Au  -wv2  +  F(l-u) 

(4-1) 

v,  =  dv Av  +  uv2  -  (F  +  k)v 

In  the  homogeneous  state  the  concentrations  do  not  change  with  respect  to  space  and  there  will 

be  no  diffusion.  In  other  words,  the  Laplacian  terms  that  represent  diffusion  go  to  0.  The 

resulting  equations  are: 

u  =  f(u,v)  =  -uv2+F(\-u) 

,  (4.2) 

v'  =  g(u,v )  =  uv 2  -  (F  +  k)v 

It  is  interesting  to  note  that  for  periodic  or  Neumann  boundary  conditions,  (4.2)  also  results  when 
considering  only  the  0th  order  Galerkin  approximation  of  (4. 1). 

Theorem  4.1:  When  using  periodic  or  Nuemann  boundary  conditions,  the  diffusionless 
equations  describing  the  reaction  kinetics  in  (4.2)  are  equivalent  to  the  ODE’s  that  result  from 
the  Galerkin  method  when  choosing  N=0. 

Proof:  When  </>m„(x,y)  =  em(mx+ny)  as  for  periodic  boundary  conditions,  or  else  when 


(j) m  n  ( x ,  y)  =  cos(;r  mx)  cos(;r  ny)  as  for  Neumann  boundary  conditions,  <f>0  0  (x,  y)  =  1 . 
solution  template  (3.46)  with  N=0  now  takes  the  form: 


0  0 

2(x,  y,t)  =  XXa«%,(^j)  =  ao,o(?M>,o  {x,y) 

m= 1  n= 1 


k,o(0 
[^0,0  (0- 


(4.3) 


The 


Now  when  we  substitute  the  right  side  of  (4.3)  into  (4.1),  we  have 

ao,o  —  duAa00  —  a00b0  0  +  A(1  —  a0  0)  =  — o0  060  0  +A(1  —  a00) 

2  2  0 4-4 ) 

Ko  =  dvAb0.o  +  ao,oK,o  —  {F  +  k)b00  =  a00b00  —(F  +  k)b00. 

Since  the  coordinate  functions  a0  0  and  b0  0  depend  only  on  time,  the  diffusion  tenns  in  (4.3) 

that  contain  the  Laplacian,  a  spatial  differential  operator,  go  to  0.  Removing  the  subscripts  we 
see  that  (4.4)  is  equivalent  to  (4.2): 


a'  =  -ab1  +F(l-a) 
b'  =  ab2  -(F  +  k)b. 


(4.5) 


42 


Thus  (4.5),  the  system  of  equations  that  arises  for  the  N=0  case,  is  equivalent  to  (4.2). 

Q.E.D. 

The  fact  that  the  equilibrium  points  we  are  looking  for  are  to  be  homogeneous  allowed  us  to 
remove  the  diffusion  term.  The  fact  that  they  are  equilibrium  points  means  that  they  will  not 
change  over  time,  and  thus  the  left  side  of  (4.5)  will  be  0.  In  order  to  find  these  homogeneous 
equilibrium  points  we  must  now  solve  the  nonlinear  system: 

-ab2  +F(l-a)  =  0 

(4.6) 

ab 2  -(F  +  k)b  =  0. 

Solving  for  a  in  the  first  equation  and  substituting  into  the  second  equation  in  (4.6),  we  see  that 
b(b2(F  +  k)-bF  +  F2 +Fk)  =  0,  (4.7) 

thus  either  b=  0,  in  which  case  <7=1,  or  else  b  will  be  the  root  of  the  polynomial 
b2  (F  +  k)  -  bF  +  F2  +  Fk  .  We  will  show  later  that  the  equilibrium  point,  (<7o,ho)=(l,0)  will 
always  be  stable,  even  when  diffusion  is  added.  The  b  coordinates  of  the  other  equilibrium 
points  are,  by  the  quadratic  formula, 
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less  than  the  coordinate  b2  and  ci\  will  always  be  greater  than  a2.  We  will  call  (ao,bo)=(  1,0)  the 
independent  equilibrium  point  and  the  other  two  dependent  equilibrium  points  since  they  depend 

on  F  and  k.  The  dependent  equilibrium  points  will  only  exist  when  k  <  4f  jl-F  .  In  that 
region  of  F-k  space,  (4.2)  has  three  equilibrium  points  (See  fig  4.1).  Elsewhere,  ( ao,bo )  is  the 
only  equilibrium  point.  This  is  illustrated  in  Figure  4.1,  which  is  also  reproduced  in  [18].  Using 
calculus  we  can  find  the  maximum  k  value  where  there  will  be  exactly  2  equilibrium  points: 


7  yfF 

k  = - F 


dk 

dF 


4  y[F 


The  critical  point  is  at  F=k=  1/1 6=0.0625  as  it 
appears  to  be  in  Figure  4.1.  The  maximum  value 
of  F  for  which  the  dependant  equilibrium  points 
can  exist  (see  Figure  4.1)  occurs  at  the  larger 
value  for  which  exactly  two  equilibrium  points 
exist  when  k  is  equal  to  0.  The  values  of  F  for 
which  there  are  exactly  two  equilibrium  points 


with  k=0  are  given  by  F  = 


F=0  03  k=0  04  time=2 

11WW 


Figure  4.2  Field  plot  and  equilibrium  points 
of  the  ODE  system  (4.2)  for  F=0.03,  k=0.04. 


Solving  for  F  we  have  F= 0  and  F=  1/4  as  indicated  in  Figure 


4.1.  Thus,  region  of  F-k  space  on  which  the  dependent  homogeneous  equilibrium  points  exist  is 
bounded  by 


0  <  F  <  1/4 
0  <  A:  <  1/16. 


(4.10) 


Figure  4.2  shows  the  equilibrium  points  and  the  vector  field  plot  of  the  differential 
equation  in  the  ab  plane  for  F=0.03  and  A=0.04.  At  these  particular  parameter  values,  there  are 
three  equilibrium  points.  Figures  4.3  and  4.4  show  the  dependence  of  the  equilibrium  points 
given  by  (4.8)  and  (4.9)  on  F  and  k  as  these  parameters  are  respectively  held  constant.  In  figure 
4.3,  the  equilibrium  points  approach  each  other  from  either  side  as  k  is  increased.  In  Figure  4.4, 
the  equilibrium  points  appear,  traverse  around  an  enclosed  path,  and  vanish  as  F  increases  from 

(1/4  -  Vl/16  -  A:)2  to  (1/4  +  V1/16-A:)2. 
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The  surface  and  the  contour  plot  in  Figure  4.5  show  how  the  equilibrium  points  change 
with  both  F  and  k.  Analyzing  Figure  4.5  we  can  visualize  several  properties  of  the  system: 


equilibrium  points  for  fixed  k,  increasing  F 


a 


Figure  4.3:  Equilibrium  points  as  k  varies  and  F  is 
held  constant. 


Figure  4.4:  Equilibrium  points  as  F  is  varied  and 
k  is  held  constant. 


Theorem  4.2:  When  there  are  exactly  two  equilibrium  points,  a=l/2. 


Figure  4.5:  Visualizing  the  dependence  of  the  equilibrium  points  on  F  and  k. 


Proof:  When  there  are  exactly  two  equilibrium  points,  the  term  under  the  radical  in  (4.8)  is  zero. 


This  is  true  when 


2 


F  .  Thus  we  see  that  the  equilibrium  point’s  a  coordinate  is  at: 
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F  +  k 

a  = - 

b 


{F  +  k) 

'  F  ' 
,2  {F  +  k), 


2  {F  +  k)2 
F 


F 


Q.E.D. 


Theorem  4.3:  lim(a  +  b)  =  1 . 

-  *->  o 


Proof: 


As  k  — »  0  we  see  from  (4.8)  that  b  — > 


F±^F2 -4F(F2) 
IF 


1  / 2 ±  Vl/4 -f7  and  from 


(4.9)  that  a  — >  .  Thus 

F  F  1/2  +  Vl/4 -F  F(l/2  +  Vl/4-F) 

limn  = - .  = - - - -  = - 

*-*  l/2±Vl/4-F  l/2±Vl/4-F  1/2  +  V1/4-F  l/4-(l/4-F) 

l-(l/2±Vl/4-F)  =  l-limh, 

*•->  0 


or  in  other  words  lim(n  +  h)  =  1 . 

/t->0 

Q.E.D. 


Now  that  we  have  a  sense  of  the  dependence  of  the  equilibria  on  the  parameters,  we  can  ask 
whether  or  not  they  are  stable.  An  equilibrium  point  of  (4.2)  is  defined  as  stable  if  small 
perturbations  from  the  equilibrium  will  eventually  settle  back  to  the  equilibrium  point.  As  is 
explained  in  depth  in  [10],  p.  178-179,  [22],  pl50,  and  in  [15],  p.  697-701,  the  stability  of  a 
nonlinear  system  can  be  detennined  from  the  eigenvalues  of  the  matrix  obtained  by  linearizing 
the  system  about  the  equilibrim  point.  Let  u=p  and  v=q  be  an  equilibrium  point  of  the  system, 
so  that: 

f{p,q)  =  g{p,q)  =  0-  (4.11) 

We  linearize  the  system  about  (p,q ).  That  is,  we  take  a  Taylor  series  expansion  of/ and  g  for  u 
and  v  very  close  to  {p,q )  and  keep  only  the  linear  terms.  We  set  a={p+s)  and  b=(q+3)  where  s 
and  d  are  small,  and  find  the  Taylor  series  expansion  about  (p,q): 

a,  =  f{a,b)  =  f(p  +  £,q  +  S)  =  f(p,q)  +  £^(p,q)  +  S^{p,q)  +  higher  order  terms 

b,  =  g{a,b)  =  g(p  +  £,  q  +  8)  =  g(p,  q)  +  £  f[  (p,  q)  +  8  %  (p,  q)  +  higher  order  terms. 


(4.12) 
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Noticing  that  a,  =st  and  bt  =St,  using  (4.11),  and  ignoring  the  higher  order  terms  of  (4.12) 

because  they  are  proportional  to  the  small  numbers  e  and  S  raised  to  the  second  degree  and 
higher,  we  rewrite  (4.12)  as: 


e,  =ei{p,q)  +  S^{p,q) 


The  solutions  to  (4.13)  depend  on  the  eigenvalues  of  the  matrix 


(4.14) 


which  we  will  call  the  stability  matrix.  The  notation  indicates  that  it  is  evaluated  at  the 
equilibrium  point,  (p,q ),  where  we  would  like  to  test  stability. 

The  solutions  to  (4.13)  will  contain  terms  with  the  eigenvalues  of  Jin  the  exponent.  Thus 
if  any  of  the  eigenvalues  of  J  has  a  positive  real  part,  the  perturbations  3  and  s  will  increase 
without  bound  and  the  system  is  unstable.  If  J’s  eigenvalues  have  negative  real  parts,  then  the 
perturbations  will  decrease  exponentially  with  time  and  the  system  is  stable. 

For  the  Gray-Scott  reaction  system  (4.2),  the  stability  matrix  is  given  by: 


J 


«oA 


-b2-F  -lab 
b 2  ab-(F  +  k)_ 


(4.15) 


From  (4.15)  we  can  see  that  the  eigenvalues  of  Jm,  and  thus  the  linear  stability  of  (4.2)  at  (p,q), 
depends  on  the  equilibrium  points  and  the  parameters  F  and  k.  From  (4.6),  (4.8)  and  (4.9),  we 
know  that  the  equilibrium  points  themselves  depend  on  the  parameter  values  F  and  k,  so  the 
stability  of  the  eigenvalues  depends  ultimately  on  F  and  k  in  our  model. 

Theorem  4.4:  The  independent  equilibrium  point,  (ao,bo)=(l,0)  is  always  stable  for  (4.2). 

Proof:  At  ( ao,bo ), 


J 


a0,bQ 


0  -(F  +  k)  J 


(4.16) 
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For  diagonal  matrices,  the  eigenvalues  are  the  diagonal  entries,  so  (4.16)  has  eigenvalues  of-i7 
and  -( F+k ).  Since  the  feed  rate  F  and  the  rate  constant  k  will  always  be  positive  number,  the 
eigenvalues  of  (4.16)  will  always  be  negative  for  any  valid  choice  of  F  and  k. 

Q.E.D. 


Theorem  4.5:  The  dependent  equilibrium  point,  ( a\,b\ )  is  always  unstable  for  all  plausible 
values  of  F  and  k. 

Proof:  At  (a\,b\). 


J 


a,  ,b, 


- b 2  -  F 
b2 


-2 


F  +  k 


F  +  k 


b- (F  +  k) 


- b 2  -  F 
b2 


-2  (F  +  k) 
F  +  k 


The  determinant  of  this  matrix  is 


Det(Jaibx  )  =  (-b2-  F)(F  +  k)  +  2(F  +  k)b2  =(F  +  k)(b2  -  F) 


(4.17) 


Since  the  entries  of  the  matrix  Ja  b  are  all  real,  the  eigenvalues  of  the  matrix  will  either  be  real 

or,  if  they  are  complex,  complex  conjugates  of  one  another.  If  we  can  show  that  the  determinant 
of  JaA  is  negative,  we  will  know  that  its  eigenvalues  are  real  and  that  one  eigenvalue  is  positive, 

while  the  other  eigenvalue  is  negative.  That  is  because  the  determinant  of  a  matrix  is  the  product 
of  its  eigenvalues.  The  factor  (F+k)  is  always  positive  because  F  and  k  are  both  positive  values. 
Our  goal,  therefore,  is  to  show  that  the  term  ( b 2  -  F)  is  always  negative  for  all  plausible  values 
ofAandk.  Note  here  that  from  (4.10)  that  0  <  F  <  1/4  and  0  <  k  <  1/16  . 

The  equilibrium  point’s  b  value,  hi,  is  given  by: 


F-Jp2 -4F(F  +  k)2 
bl  ~  2 (F  +  k)  ' 

The  term  inside  the  radical  must  be  between  0  and  1 ; 

0  <F2 -4 F(F  +  k)2  <1 

since  it  must  be  positive  for  hi  to  be  real  and  since  F2<F<1/4<1.  It  then  follows  that: 


F2  -  4 (F  +  k )2  <  F2  -  4 F(F  +  k )2  <  ^F2  -4F(F  +  k)2  . 
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The  left  inequality  is  true  because  F  <  1  /  4  <  1  and  the  right  inequality  is  true  because  the  term 
inside  the  radical  is  less  than  one.  We  rearrange  the  resulting  inequality  to  get  b  on  the  left  side: 


F2  -  ^F2  -4  F(F  +  k)2  <  4 (F  +  kf 


F2  -Jf2  -4F(F  +  k)2 

b  = - - - - - —  <2  (F  +  k) 

2  (F  +  k) 


Now  we  manipulate  this  expression  to  get  it  into  a  more  convenient  form: 


bF 


(F  +  k) 


<2  F 


bF-F(F  +  k) 
(F  +  k) 


-F<  0. 


(4.18) 


2 

From  the  original  polynomial  (4.7),  we  can  get  an  expression  for  F  in  terms  of  b\ 


bF-F(F  +  k) 
(F  +  k) 


(4.19) 


Now  using  (4.18)  and  (4.19)  we  can  show: 


y  F  _bF  ~F(F  +  k) 
(F  +  k) 


F<  0, 


which  means  that  the  determinant  of  the  matrix  J  ^  is  negative.  Therefore  both  its  eigenvalues 

are  real  and  one  eigenvalue  is  positive,  while  the  other  is  negative.  The  equilibrium  point  (a\,b\) 
will  always  be  unstable  for  valid  values  of  F  and  k. 


Q.E.D. 
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Determining  stability  at  the  dependent  equilibrium  points  {a.2,b 2)  is  slightly  more  complicated 
than  for  {a\,b\)  and  (ao,bo)  since  it  has  complex  eigenvalues.  Figure  4.6  is  a  plot  of  the 


Eigenvalues  of  the  stability  matrix  of  (avb7) 


Fig.  4.6  Eigenvalues  of  the  stability  matrix  for  fixed  k  as  F  varies 


E>o*nvafaM  off*  lUDttj  m«*»  of 


eigenvalues  of  (4.15)  in  the  complex  plane  evaluated  at  the  equilibrium  points  as  F  varies 
between  0.01  and  0.17  and  k  is  held  at  0.4.  The  eigenvalues  of  the  dependent  equilibrium  point 
(02,62)  are  the  complex  conjugates  of  one  another,  and  they  cross  into  the  left  half  plane  with 
nonzero  imaginary  components.  When  this  occurs,  a  family  of  periodic  solutions  appear 
surrounding  the  equilibrium  point.  This  phenomenon  is  known  as  Hopf  bifurcation  and  is 
discussed  in  detail  in  [15],  p  706-719.  The  dotted  line  in  Figure  4.7  indicates  where  in  F-k  space 
these  Hopf  bifurcation  points  will 
occur.  Notice  the  behavior  of  the 
eigenvalues  in  Figure  4.6  in  the 
different  portions  of  the  vertical  line 
at  6=0.04.  When  we  add  diffusion 
we  will  search  for  patterns  around 
this  Hopf  bifurcation  line  in 
parameter  space  reasoning  that  these 
periodic  solutions  in  the  N=0  case 


The  dotted  line  in  Figure  4.7  indicates  where  in  F-k  space 
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might  give  rise  to  nonhomogeneous  steady-state  solutions  and  periodic  orbits  when  diffusion  is 
added.  Figure  4.7  is  a  reproduction  of  a  similar  diagram  found  in  [18]. 

From  Figure  4.6  we  can  also  see  that,  as  we  established  earlier,  ( ao,bo )  is  always  stable, 
and  (ci\,b\)  is  always  unstable.  The  fact  that  (a\,b\)  has  one  negative  and  one  positive  eigenvalue 
also  means  that  it  is  a  saddle  point.  Figure  4.8  shows  the  phase  plane  of  (4.2)  for  F=0.011, 
&=0.04.  These  parameter  values  lie  between  the  dotted  line  and  the  lower  solid  line  in  Figure 
4.7.  At  these  values  we  will  therefore  have  three  equilibrium  points,  two  of  which  will  be 
unstable.  The  small  cross  in  Figure  4.8  indicates  the  location  of  (a\,b\),  and  the  circles  represent 
the  other  equilibria.  Since  (a\,b\)  is  a  saddle  point,  it  must  have  an  unstable  and  a  stable 
manifold.  That  is,  there  must  be  a  curve  along  which  solutions  are  attracted  and  a  curve  along 
which  solutions  are  repelled.  These  are  shown  in  the  figure,  and  they  help  us  understand  why 
solutions  arising  from  initial  values  very  close  to  the  stable  equilibrium  point  ( ao,bo )  sometimes 
follow  roundabout  trajectories  before  they  finally  settle  into  (ag,bo). 


Stability  Diagram  for  a  parameter  value  with  three  equilibria,  two  of  which  are  unstable 
F=0.011  k=0.04  time=1 00 


Fig  4.8  Stable  and  unstable  manifolds  of  the  saddle  point  (aub i) 

Now  that  we  understand  the  behavior  of  the  system  without  diffusion,  we  will  try  to  get  a 
sense  of  what  happens  once  diffusion  is  added.  We  will  first  derive  the  conditions  for  linear 
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stability  in  the  system  once  the  diffusion  term  is  included  in  the  analysis.  Consider  the  general 
two  species  reaction-diffusion  equation: 
u,  =  d„Au  +  f  (u,v), 

‘  (4.20) 

vt  =  dvAv  +  g(u,v). 


We  apply  the  standard  stability  analysis  of  partial  differential  equations  by  analyzing  the 
behavior  of  solutions  of  (4.20)  in  a  small  neighborhood  of  ( uo,vo ).  This  procedure  is  similar  to 
the  stability  analysis  outlined  above.  To  that  end,  we  make  the  substitutions 

U  =u0  +  £U, 

V  =  Va+£V. 


Linearizing/ and  g  near  the  equilibrium  point  (uo,v0)  by  using  a  Taylor  series  expansion  about 
(uo,vo),  we  find  that  for  small  s: 


f(U,  V )  =  f(u0  +£,v0+s)  =  f(u0  ,v0)  +  £ 


g(U  ,V)  =  g(u0  +£,v0  +£)  =  g(u0,V0)  +  £ 


ydU 

3gj 

du 


d[_ 

1  ("O  >V0  )  ^  Qy 

,  dg  , 

u  ~dv 


( «o .v0 ) 


+  0(e  )  +  ..., 


+  0(£2)  +  ...  . 


But  since  we  are  concerned  only  with  very  small  s  the  0(s")  terms  can  be  ignored.  Also  since 
(«o,vo)  is  an  equilibrium  point,  the  constant  terms,  /(w^'’o)  =  g(uo,  Vo)  are  zero.  We  are  thus  left 
with: 


f(U,V)*£ 

g(U,V)~£ 


Kdu 

(  dg 


u  + 


u  + 


V 

dv 

dg 

dv 


l(«0 -vo) 


l(«0.v0) 


=  £j 

\  u0,v0 


Vv2 


where  J  is  the  stability  matrix 

u0 ,v0  j 


~df 

df  , 

a  b 

du  l(“»'Vo> 

dv  <Uo-Vo> 

c  d 

dg  | 

dg  , 

_du  l(“0'Vo) 

dv  (l,°’Vo). 

(4.21) 


of  the  functions / and  g  evaluated  at  ( u0,vq ).  We  also  note  that 


U,  =£Ut 

V,  =£V, 


A  U  =£Au 
AV  =  £  Av  ■ 


and 
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Substituting  U  and  V  into  (4.20),  dividing  out  by  s,  and  ignoring  the  remaining  terms  that  depend 
on  s,  we  see  that  near  equilibrium  a  and  v  satisfy  the  following  linear  partial  differential 
equation: 

u,  =  dAu  +  au  +  bv, 

'  "  (4.22) 

v,  =dvAv  +  cu  +  dv. 


Now  let  ( uo,vo )  be  a  linearly  stable  equilibrium  point  for  the  diffusionless  system 
u,  =f(u,v), 

V,  =g(u,v). 

Note  that  the  linearization  of  (4.23)  about  (uo,vq)  is 


(4.23) 


\vtJ 


=  J„ 


U  o,v0 


\vj 


The  assumption  that  (mo,vo)  is  linearly  stable  is  equivalent  to  assuming  that  the  real  parts  of  the 
eigenvalues  of  J  v  are  always  negative.  Let  k\  and  kn  be  the  eigenvalues  of  J u  v  .  The  fact  that 


they  are  always  negative  means  that: 

a  +  d  =  Tr(J)  =  A,  +  A2  <  0 ,  (4.24) 

and 


ad  -  be  =  Det(J)  =  TjX,  >  0 .  (4.25) 

We  seek  solutions  to  (4.14)  of  the  form 


u(x,y,t)  =  AeA,(ptj, 
v(x,y,t)  =  Be/'(piJ, 


(4.26) 


where  <pj(x,y )  is  an  eigenfunction  of  the  laplacian  operator  A,  i.e.,  A  <Pi/x,y)=  Sy  <Pj(x,y)  where  by 
is  a  negative  scalar  depending  on  the  particular  basis  function  used.  Substituting  (4.26)  into 
(4.20)  yields 

=  SijAdueA'<pij  +  aAe'J<Pj  +  bBeXt  (ptj 
XBeh (ptj  -8yBdveM(py  +cAek(py  +  dBe/J (pir 


After  dividing  by  ek<Pj  we  have: 
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~A 

a  +  Sydu 

b 

~A 

B 

c 

d  +  $  yd  v _ 

B 

So  for  solutions  described  by  (4.26)  to  exist,  X  must  be  an  eigenvalue  of  the  matrix 


Dij  = 


a  +  djjdu  b 


d  +  Sjdv 


(4.27) 


In  order  for  (4.26)  to  be  a  linearly  stable  solution  of  (4.20),  the  eigenvalues Xa  and/lfe  of  (4.27) 
must  all  have  negative  real  parts.  Using  the  properties  of  determinants  and  trace,  this  means  that 
Sy2(dudv )  +  Sy  (ddu  +  adv )  +  {ad -be)  =  Det(Dy  )  >  0 ,  (4.28) 

and 


a  +  d  +  Sy  (du  +dv)  -  Tr{Dy )  -  Aa  +  Ab  <0 .  (4.29) 

These  are  necessary  and  sufficient  conditions  for  stability,  whether  or  not  the  eigenvalues  are 
complex. 

One  phenomenon  that  is  of  interest  in  pattern  formation  in  reaction-diffusion  systems  is 
that  of  diffusion  driven  instability.  A  homogeneous  equilibrium  point  that  exhibits  diffusion 
driven  instability  will  be  stable  for  the  diffusionless  system  (4.23),  but  will  become  unstable  once 
diffusion  is  accounted  for.  Diffusion  is  usually  a  process  that  contributes  to  stability  in  a  system 
by  decreasing  concentrations  where  they  are  unusually  high,  and  increasing  concentrations  where 
they  are  lower  than  the  surrounding  regions.  In  our  reaction-diffusion  system,  we  have  found  a 
region  of  parameter  values  for  which  diffusion  breaks  the  linear  stability  of  a  system.  We  have 
found  evidence  of  stable  periodic  and  nonhomogeneous  steady-state  solutions  arising  in  this 
parameter  region  as  the  systems  proceed  away  from  the  unstable  homogeneous  equilibrium. 
Murray  and  others  ([15], [17])  have  shown  in  other  reaction-diffusion  systems  that  pattern 
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formation  in  the  form  of  periodic  solutions  and  nonhomogeneous  steady-state  solutions  can 
occur  for  parameter  values  regions  where  diffusion  driven  instability  has  been  shown  to  exist. 
Perhaps  the  instability  brought  on  by  the  normally  restorative  process  of  diffusion  combined  with 
the  fact  that  the  system  has  several  stable  homogeneous  solutions  without  diffusion  is  enough  to 
make  the  system  unstable  enough  to  drive  solutions  away  from  the  homogeneous  state,  while 
remaining  stable  enough  to  find  the  nonhomogeneous  steady-state  and  periodic  solutions  in 
between  homogeneous  equilibria  that  we  call  patterns.  The  following  result  is  a  necessary 
condition  for  diffusion  driven  instability. 

Theorem  4.6:  In  order  for  diffusion  driven  instability  to  occur  in  a  two  species  reaction- 
diffusion  equation,  the  diffusion  constants  must  be  unequal  In  other  words,  in  order  for  a 
homogeneous  equilibrium  solution  ( ui,vi )  to  be  linearly  stable  for  the  diffusionless  system  (4.23), 
and  linearly  unstable  for  the  full  system  (4.24),  du  +  dv . 

Proof:  In  order  for  the  homogeneous  equilibrium  solution  (u\,v\)  to  be  linearly  stable  for  (4.23), 
conditions  (4.24)  and  (4.25)  must  hold,  i.e.  a+d< 0  and  ad-bc> 0.  If  we  let  the  diffusion  constants 
equal  one  another,  du=dv=duv,  then  the  left  side  of  condition  (4.28)  becomes 
Sy 2dJ  +  bljdufa  +  d)  +  (ad -be)  . 

The  first  term  in  this  expression  will  always  be  positive.  The  second  term  is  also  positive  since 
Sj  is  a  negative  number,  the  diffusion  constant  is  a  positive  number,  and  (a+d)  is  negative.  The 

last  tenn  is  also  positive,  so  (4.28)  always  holds. 

With  equal  diffusion  constants,  the  left  side  of  (4.29)  becomes: 

a  +  d  +  2Sijduv , 
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which  by  similar  reasoning  is  always  negative  causing  (4.29)  to  hold  for  all  parameters.  Thus 
always  satisfies  sufficient  conditions  for  stability  in  the  full  system  (4.22)  and  will  always 
be  stable.  There  can  be  no  diffusion  driven  instability. 

Q.E.D. 

We  will  now  see  how  this  condition  can  be  applied  to  the  Gray-Scott  model  to  find  regions  of 
diffusion  driven  instability.  The  first  thing  we  would  like  to  show  is  that  the  homogeneous 
equilibrium  point  ( ao,bo )  remains  linearly  stable  for  any  choice  of  F,  k  and  diffusion  constant  du 
and  dv. 


Theorem  4.7:  The  homogeneous  equilibrium  solution  to  the  Gray-Scott  equations  (4.1)  at 
(wo,vo)=  (1,0)  is  linearly  stable  for  all  feasible  values  ofF,  k,  du  and  dv. 

Proof:  First  of  all,  we  know  that  («0,vo)=(l,0)  is  a  homogeneous  solution  for  (4.1),  because  the 
point  (ao,bo)=(l,0)  is  an  equilibrium  point  of  the  diffusionless  system.  In  fact,  Theorem  4.4 
states  that  ( ao,bo )  is  linearly  stable  for  all  F  and  k.  For  the  Gray-Scott  reaction-diffusion 
equation,  the  matrix  /from  (4.21)  is 


a  b 

~-v2  -  F 

-2  uv 

c  d 

v2 

2  uv  —  (F  +  k) 

and  the  matrix  Dy  is  given  by 


Dij  = 


a  +  Sijdu  b 
c  d  +  8 jd v 


vz-F  +  Sjdl, 


-2  MV 

2  mv  -(F  +  k)  +  8ijdv 


.  (4.30) 


Substituting  the  values  (wo,vo)=(l,0)  into  (4.30)  we  have 


a  +  Sjdu  b 

-F  +  8 yd  u 

0 

c  d  +  8iJdv 

0 

-(F  +  k)  +  8ydv 

(4.31) 


The  left  side  of  condition  (4.28)  becomes: 


8„2dady  -  S0 ( Fdv  +  (F  +  k)du )  +  F(F  +  k ) . 
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Since  Sy  is  always  negative  and  the  reaction  and  diffusion  parameters  are  always  positive,  this 

expression  is  always  positive  and  (4.28)  always  holds. 

For  (cio, bo),  the  left  side  of  condition  (4.29)  becomes: 

-(2  F  +  k)  +  Sy(du+dv). 

This  expression  is  always  negative,  and  thus  condition  (4.29)  holds.  Since  both  conditions  hold, 
the  equilibrium  point  (do, bo)  is  stable  for  all  valid  reaction  and  diffusion  parameter  values. 

Q.E.D. 


For  (u2,V2)=(d2,b2),  the  other  equilibrium  point  that  can  be  linearly  stable  for  the 
diffusionless  system,  the  conditions  for  stability  in  the  full  system,  (4.28)  and  (4.29),  may  not 
always  be  satisfied.  In  order  to  examine  the  stability  of  the  full  system  with  diffusion,  we  will 
look  at  the  real  parts  of  the  eigenvalues  of  the  stability  matrix.  In  the  Hopf  bifurcation  region, 
when  (02,62)  first  becomes  stable,  the  eigenvalues  will  always  be  complex  conjugates  of  one 
another.  Thus  the  real  part  of  one  eigenvalue  will  be  equal  to  the  other.  If  this  value  is  negative 
the  system  is  stable,  and  if  it  is  positive  the  system  is  unstable.  Figure  4.9  is  a  plot  of  the  real 
part  of  one  eigenvalue  of  the  matrix  Dy  (4.30)  for  F=0.0154,  k=0.04,  dv=0. 00001,  and  du  values 
of  0.0001  and  0.00002  for  the  homogeneous  equilibrium  point  (<22,62)  versus  Sy  . 

Figure  4.9  Real  parts  of  the  eigenvalues  of  matrix  (4.30)  as  is  increased.  The  plot  on  the  left  has 
</„=0. 00002  and  the  plot  on  the  right  has  </„=0.0001.  All  other  parameters  are  equal  at  .F=0.0154,  k=0.04,  and 


</„=0.00001. 
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The  plot  on  the  left  shows  that  the  real  part  of  the  greatest  eigenvalue  of  the  stability  matrix 
never  becomes  positive,  and  thus  diffusion  driven  instability  does  not  occur.  On  the  right,  we 
see  that  the  real  parts  of  the  initially  stable  point  start  negative,  but  suddenly  one  real  part  breaks 
up  and  becomes  positive.  At  this  point  the  eigenvalues  are  no  longer  complex  and  the  real  part 
of  the  other  eigenvalue,  which  remains  negative,  is  not  shown.  The  system  becomes  linearly 
unstable  and  diffusion  driven  instability  exists.  We  will  look  at  patterns  for  both  of  these 
diffusion  parameters,  and  see  what  effect  this  type  of  instability  seems  to  have  on  the  system. 
This  analysis  seems  to  indicate  at  least  initially  that  the  magnitude  of  the  difference  between  the 
diffusion  parameters  of  the  morphogens  is  a  primary  factor  in  diffusion  driven  instability. 

Murray  and  others  have  examined  diffusion  driven  instability  searching  for  Turing 
patterns.  Turing’s  original  idea  was  that  patterns  could  be  driven  by  the  diffusion  terms  and  not 
the  reaction  terms  [15],  [17],  [25].  We  have  found  some  indications  of  pattern  formation  inside  an 
envelope  of  parameter  values  where  diffusion  driven  instability  can  exist  for  certain  modes  and 


diffusion  coefficients. 
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Results 

This  section  displays  some  of  the  patterns  that  we  observed  for  various  parameter  values.  The 
values  of  the  reaction  parameters  that  we  used,  F=0.0154  and  A:=0.04,  fall  above  the  Hopf 
bifurcation  curve  of  Figure  4.7.  With  N=0,  therefore  (<22,62)  ~  (0.2748,  0.2016)  would  be  a 
stable,  spiral  equilibrium  point.  Figure  5.1  shows  snapshots  of  our  simulated  solution  of  u  at 
these  parameter  values  with  N=2.  Recall  that  our  diffusion  constants  are  du=2  x  1 0°  and  dv=  10°. 
Our  initial  condition  is  a  depression  centered  at  (x,y)=(-0.65,  0.65).  By  the  time  t  reaches  5801, 
we  can  see  the  form  of  a  morphogen  wave  which  becomes  fully  established  by  time  9801. 

Figure  5.1  Simulated  solution  for  u  of  the  Gray-Scott  model  with  N=2. 


Gray-Scott  model,  n=2  du=.00002,  dv=.00001 ,  k=.04,  F=.0154,  t=1  Gray-Scott  model,  n=2  du=.00002,  dv=. 00001 ,  k=.04,  F=0154,  t=1001 


Gray-Scott  model,  n=2  du=.00002,  dv=,00001 ,  k=.04,  F=. 0154,  t=1801 


Gray-Scott  model,  n=2  du=.00002,  dv=,00001 ,  k=.04,  F=.0154,  t=3801 


59 


Gray-Scott  model,  n=2  du=  00002,  dv=.00001 .  k=.04,  F=.0154,  t=5801 


1 


Gray-Scott  model,  n=2  du=.00002,  dv=. 00001 .  k=.04,  F=.0154,  t=9801 


In  Figure  5.2  we  have  plots  of  a  few  of  the  coordinate  functions  of  time  from  the  solution 
template  (3.59),  the  ald(t)  and  bj/t)  versus  time  for  the  solution  from  Figure  5.1.  Notice  that  the 
imaginary  parts  for  ao,o(t)  and  ho,o(0  are  essentially  0.  Near  time  6000  we  see  that  the  coordinate 
functions  for  the  0,0  modes,  which  we  will  call  the  principal  coordinate  functions,  seem  to 


Figure  5.2  -  Several  coordinate  functions  of  time  for  F=0.0154,  &=0.04,  N=2.  Real  parts  are  on  the  right  and 
imaginary  parts  are  on  the  left.  The  captions  on  the  individual  pictures  refer  to  the  index  of  the  a  vector  in 
the  MATLAB  program  GAL  PEARSON2D  (see  appendix). 
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achieve  a  steady  state.  The  values  they  obtain  are  not  far  from  the  spiral  homogeneous 
equilibrium  of  the  N=0  problem  (<22, b?)  ~  (0.2748,  0.2016).  The  coordinate  function  ai, i{t) 
begins  to  oscillate  steadily  around  this  same  time.  This  time  also  corresponds  to  the  time  at 
which  the  wave  begins  to  appear  in  Figure  5.1. 

Figure  5.3  shows  what  happens  when  we  increase  N.  We  plot  u  for  the  same  parameter 
values,  except  that  N=4  for  times  up  to  1000.  This  time  we  have  placed  the  initial  perturbation 
in  the  center,  but  this  does  not  appreciably  affect  the  results.  The  perturbation  spreads  to  the 
boundaries  and  from  the  site  of  the  perturbation  a  wave  seems  to  be  emanating  by  time  1000. 
Figure  5.4  shows  the  contour  plots  of  u  with  N=4,  this  time  with  the  initial  perturbation  off 
center  again.  As  time  increases  past  1000,  the  circular  waves  emanating  from  the  site  of  the 
initial  perturbation  seem  to  take  on  transverse  wavelike  motion,  similar  to  those  we  saw  in  Figure 
5.1  but  not  as  clearly  defined. 

As  N  is  increased,  we  are  able  to  use  more  modes  to  approximate  the  solution  and  thus 
we  are  able  to  capture  more  detail  and  complexity.  If  waves  form  in  the  true  solution  for  these 
parameter  values,  they  may  not  look  like  those  we  obtained  for  N=  2  and  N= 4.  However,  if  the 
true  solutions  do  indeed  have  waves,  these  low  level  approximations  could  prove  useful  in 
predicting  the  formation  of  waves  for  a  particular  parameter  set.  Comparing  the  solutions  near 


time  1000  for  n= 2  and  n= 4,  there  is  little  resemblance.  We  cannot  be  confident  that  our 
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simulated  solution  is  close  to  the  true  solution  until  we  can  increase  the  number  of  modes 
without  any  significant  change  in  the  simulated  solution.  So  far  N=4  is  the  highest  number  of 
modes  we  have  been  able  to  use  with  MATLAB.  If  we  were  to  use  Neumann  boundary 
conditions  and  a  smaller  region,  we  might  be  able  to  try  more  modes. 
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Figure  5.3  Surface  Plots  of  iV=4  with  FMJ.0154,  A=0.04,  for  times  up  to  1000. 


Gray  Scott  model,  n=4  du=.00002,  dv=. 00001 ,  k=.04,  F=.0154,  initial  condition 


Gray  Scott  model,  n=4  du=. 00002,  dv=. 00001 ,  k=.04,  F=.0154,  t=60 


Gray  Scott  model,  n=4  du=. 00002,  dv=. 00001 ,  k=.04,  F=.0154,  t=1 40 


Gray  Scott  model,  n=4  du=. 00002,  dv=.00001 ,  k=.04,  F=.0154,  t=220 


Gray  Scott  model,  n=4  du=. 00002,  dv=.00001 ,  k=.04,  F=. 0154,  t=300 


Gray  Scott  model,  n=4  du=.00002,  dv=.  00001 ,  k=.04,  F=  0154,  t=380 


Gray  Scott  model,  n=4  du=. 00002,  dv=. 00001 ,  k=.04,  F=.0154,  t=420 


Gray  Scott  model,  n=4  du=. 00002,  dv=. 00001 ,  k=.04,  F=.0154,  t=1000 
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Figure  3.8  u  for  4  basis  functions  at  various  times,  f^O.0154,  k= 0.04 
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All  of  the  patterns  above  have  been  periodic,  oscillating  with  time  in  a  predictable 


manner.  Motivated  by  our  analysis  of  diffusion  driven  instability,  we  alter  the  diffusion 


Figure  5.5-  Gray-Scott  Model  with  new  diffusion  rates,  d,  =0.0001,  </,,=0.00001,  AM,  /Mi.0154, 
k= 0.04,  time  up  to  1000.  The  contour  maps  all  represent  the  concentration  of  u  where  red  is  most 
concentrated  and  blue  is  least  concentrated  except  for  the  final  picture  which  shows  the 
concentrations  of  v  near  time  1000  right  beside  its  corresponding  concentration  u. 
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constants  to  obtain  what  seems  to  be  a  nonhomogeneous  steady  state  solution.  Figure  5.5 
shows  our  simulated  solution  for  N=4,  F=0.0154,  and  AH).  04.  This  time  our  diffusion 
constants  are  du  =  10"4  and  dv  =  1 0"\  Now  the  diffusion  constants  differ  by  a  multiple  of 
10  rather  than  2  as  they  did  in  Figures  5. 1-5.4.  We  start  with  a  perturbation  in  the  center 
this  time.  The  initial  condition  on  v  is  small  positive  perturbation  from  zero  centered  at 
the  origin  as  well.  As  the  simulation  proceeds,  we  see  that  the  concentration  of  u  initially 
falls  off  from  around  the  perturbation.  Near  the  center,  u  forms  a  region  of  high 
concentration  and  then  seems  to  break  into  four  regions  of  increased  concentration  that 
spread  towards  the  boundaries.  This  process  of  forming  a  concentrated  region  near  the 
center  and  then  breaking  off  repeats  itself  until  the  solution  begins  to  settle  into  a 
checkerboard  pattern  of  alternating  high  and  low  concentrations  which  then  become 
stationary  waves.  Notice  that  the  concentration  of  morphogens  u  and  v  are  nearly 
complimentary  at  the  final  time.  Gierer  and  Meinhardt  showed  in  their  activator/inhibitor 
reaction-diffusion  system  that  the  diffusion  constants  must  differ  significantly  in  order  for 
patterns  to  occur  [10]. 

We  have  shown  numerical  evidence  to  support  the  notion  that  reaction-diffusion 
equations  have  periodic  solutions  and  nonhomogeneous  steady  state  solutions.  We  have 
only  used  a  few  modes,  and  yet  we  have  generated  some  complicated  patterns  and 
structures  for  the  Gray-Scott  model.  Other  researchers  have  been  able  to  generate  and 
classify  many  of  the  patterns  that  can  be  generated  by  various  reaction-diffusion  systems, 
([5],  [9],[1 1],[12],[14],[15],[17],[18],[19],[20],[24],[26],[28]).  These  systems,  which  are 
derived  from  the  most  basic  physical  and  chemical  laws,  seem  to  take  on  a  life  of  their 
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own  given  the  right  conditions.  As  research  in  this  area  proceeds,  it  will  be  interesting  to 
see  how  well  we  can  predict  and  control  these  self-organizing  systems. 
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Galerkin  Method 

BASIS 

function  z=basis (n, x, bas) 

%this  function  sets  up  a  basis  function  for  the  cosine  basis,  the 
chebyshev  basis,  or  the  sin  basis, 

%the  arguments  are  n-  the  integer  corresponding  to  the  basis  function  n 
is  any  integer  from  zero  up. 

%x  is  the  independant  variable,  in  symbolic  or  numeric  form. 

%basis  is  one  of  three  arguments: 

%1  gives  the  basis  cos (n*pi*x) 

%2  gives  the  basis  sin(n*pi*x) 

%3  gives  the  chebyshev  polynomial  basis  T(n) 

if  bas==l 

%cos  basis 
z=cos (n*pi*x) ; 
elseif  bas==2 

%sin  basis 
z=sin (n*pi*x) ; 
elseif  bas==3 

%chebyshev  polynomial  basis 

%it  works  symbolically,  or  numerically.  basis 

%symbolic  clause 
if  isa (x, ' sym ' ) ==1 
n=n+l ; 

%T  (0) 
if  n==l 

z=l  ; 

%T  (1) 
elseif  n==2 
z=[l,0] ; 

%T(n)  n>=3 
else 

%first  it  does  it  with  the  vectors  and  then  converts 
%to  syms 

%polynomials  represented  as  vectors 
znl- [ 1 ] ; 

z= [ 1 , 0 ]  ; 
for  j=3:n 

zn2=znl ; 
znl=z ; 

z=conv ([2,0], znl); 

z=z- [zeros (1, length ( z ) -length ( zn2 ) ) , zn2 ] ; 

end 

%now  converting  the  polynomial  back  to  syms 
m=length  (  z ) ; 

for  k=m-l : -1 : 0 

zx (m-k) = [xA ( k) ] ; 

end 

z=zx . *  z ; 
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%this  last  line  converts  a  vector  of  syms  to  the  actual 

polynomial 

z=sum ( z ) ; 
end 
return 


else 

%this  loop  does  the  chebyshev  polynomial  at  a  point  x 
n=n+l ; 
if  n==l 
z=l  ; 

elseif  n==2 

z=x; 

else 

zz (1, : ) =ones (size (x) ) ; 
z  z ( 2 ,  : ) =x ; 
for  j=3:n 

zz  ( j , : ) =2*x. *zz (j-1, : ) -zz ( j-2,  :); 

end 

z=zz  (n, : ) ; 
end 

end 


else 

'Error-not  a  valid  basis  number.  Try  help  on  BASIS  to  see  what 
numbers  are  valid' 
return 

end 

GAL  ALLENCAHNQUICK 

function  [u, x, t ] =gal_allencahnquick (nn, eps2 , tvec, xvec) 

% [u, x, t] =gal_allencahnquick (4, .01,100,  [-1 : .01:1] ) ; 

%gal_pdedisplay (u,  x,  t)  ; 

%does  the  alien  cahn  equation  ut=uxx+u-uA3  on  -1..1 
%with  boundary  conditions  u(-l,t)=-l,  u(l,t)=l, 

%and  initial  condition  u (x,  0) = . 53*x+ . 47 *sin (-1 . 5*pi*x) 

%setting  the  global  variables  for  the  nonlinearities  and  nonhomogeneous  parts 

for  use  in  the  deqn 

global  nonhomogeneous 

global  nonlinearl 

global  nonlinear2 

global  nonlinear3 

global  eps 

eps=eps2 ; 

tic, 

aa=min (xvec) ; 
bb=max (xvec) ; 
for  j  =1 : nn 

nonhomogeneous ( j ) =quadl ( @nonhomog, aa, bb,  [  ] ,  [  ]  ,  j  )  ; 
initialcond ( j ) =quadl (@init, aa, bb,  [ ] ,  []  ,  j  )  ; 
for  k=l:nn 

nonlinearl ( j , k) =quadl (Snonlinl, aa, bb,  [],[],  j  ,  k) ; 
for  m=l : nn 

nonlinear 2 ( j , k, m) =quadl (@nonlin2, aa, bb, [],[], j , k, m) ; 
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for  q=l : nn 

nonlinear 3 ( j , k,m, q) =quadl ( @nonlin3, aa, bb, [ ] , [  ]  ,  j  ,  k,m,  q)  ; 
end 

end 

end 

end, 

toe 

'global  variables  done' 
tic, 

[t, a] =ode45 ( ' gal_allencahnquickdeqn ' , tvec, initialcond) ; 

'odes  solved', 
toe 

tic, 
x=xvec; 
m=length (t) ; 
for  j  =1 : nn 

sinvec (j, : )=sin(j*pi*x) ; 

end 

size  ( sinvec) 
size (x) 
size (a) 
size (a (1,  : )  ) 
for  j  =1 :m 

u ( j  ,  : )  =x+a ( j , : ) * sinvec; 
end, 
toe 


%functions  for  the  global  variables  for  the  nonlinearities  and  nonhomogenous 
parts 

function  y=nonhomog  (x, j ) 
y= (x-x. A3) ,*sin(j*pi*x) ; 

function  y=init(x,j) 

%y= (.53*x+.47*sin(-1.5*x*pi)-x) ,*sin(j*pi*x) ; 
y=(.l*x+.9*sin(-1.5*x*pi)-x) ,*sin(j*pi*x) ; 
function  y=nonlinl (x, j , k) 
y=x ,A2.*sin(j*pi*x) ,*sin(k*pi*x); 

function  y=nonlin2 (x, j , k, m) 

y=x ,*sin(j*pi*x) ,*sin(k*pi*x) ,*sin(m*pi*x); 
function  y=nonlin3 (x, j , k, m, q) 

y=sin(j*pi*x) ,*sin(k*pi*x) ,*sin(m*pi*x) ,*sin(q*pi*x); 


GAL  ALLENCAHNQUICKDEQN 

function  dadt=gal  allencahnquickdeqn (t, a) 

%dif ferential  equation  called  in  gal  allencahnquick 

global  nonhomogeneous 

global  nonlinearl 

global  nonlinear2 

global  nonlinear3 

global  eps 

nn=length (nonhomogeneous) ; 
for  j  =1 : nn 
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nonlin2vec ( j ) =0; 
nonlin3vec ( j ) =0; 
for  k=l : nn 

for  m=l : nn 

nonlin2vec ( j ) =nonlin2vec (j)+a(k)*a(m) * nonlinear 2 ( j , k, m) ; 
for  q=l : nn 

nonlin3vec ( j ) =nonlin3vec ( j ) +a (k) *a (m) *a (q) *nonlinear3 ( j , k,m, q 

)  ; 

end 

end 

end 

end 

dadt=a-eps*pi A2* ( ( 1 : nn)  '  .  A2 )  . * a+nonhomogeneous ' -3*nonlinearl*a- 
3*nonlin2vec ' -nonlin3vec ' ; 

GAL  BURGERQUICK 

function  [t,x,u]=gal  burgerquick (NN, eps, tvec, xvec) 

%  this  function  finds  a  galerkin  approximation  to  burgers'  equation 
%  ut=eps*uxx-u*ux 

%for  dirichlet  boundary  conditions  u (0, t) =u (1, t) =0  for  x=[0,l] 

%eps  is  the  inverse  of  the  reynold's  number. 

%sin(j*pi*x)  is  the  basis  function  used 
%the  ODE  solver  used  is  ODE45 

O, 

o 

%inputs:  NN--  Number  of  basis  functions  used 

%  eps--  diffusion  parameter,  inverse  of  reynold's  number 

%  tvec--  time  of  simulation-  can  be  either  a  final  time  or  a 

vector  of  times 

%  at  which  the  solution  is  evaluated 

o, 

o 

o, 

o 

%outputs:  t--  times  at  which  solution  is  simulated,  equal  to  tvec  or 

else  output  of  ode45 

%  x —  points  along  x  axis  where  x  is  evaluated,  equal  to  xvec 

%  u —  output,  values  of  u  at  points  along  x  and  t.  u  has 

length (t)  rows  and  length (x)  columns 

O, 

o 

%Syntax : 

%  [t,  x,  u] =gal_burgerquick (32,  .01,1,  [0:  .01:1] ) 

o, 

o 

%See  also:  GAL_BURGERQUICKDEQN, GAL_BURGERDISPLAY,  GAL_BURGER, 
GAL_BURGERDEQN 

tic, 

%find  initial  values  of  a's 
initf cn (1,2) 
for  j  =1 : NN 

init ( j ) =quadl (@initfcn, 1/4, 1/2,  [],  [  ]  ,  j  )  ; 

end 

'initial  condition  finished', 
toe 

%create  P  matrix  where 

P (i,  j  ,  m)  =  (sin (j  *pi*x) *sin (k*pi*x) *pi*cos  (m*pi*x)  ) 

%uses  a  trick 


global  P2  eps2 
tic, 

for  j  =1 : NN 

for  k=l : NN 

for  m=l : NN 
if  j+m==k 
P ( j , k, m) =-pi*k/ 4 ; 
elseif  abs(j-m)==k 
P ( j , k, m) =pi*k/ 4 ; 
end 

end 

end 

end 

' P  matrix  finished', 

P2=P; 

eps2=eps 

toe 

tic, 

%do  ODE  routine 

[t,  a] =ode45 ( ' gal_burgerquickdeqn ' , tvec, init) ; 

' odes  solved ' , 

toe 

tic, 

%reconstruct  solution 
x=xvec; 

for  j=l : length (t) 

u ( j , : ) =zeros (size (x) ) ; 
for  k=l : NN 

u ( j , : ) =u ( j , : ) +a ( j , k) *sin (k*pi*x) ; 

end 

end, 

'solution  constructed', 
toe 

function  z=initfcn (y, m) 
z=sin(m*pi*y) .* (1-cos (8*pi*y) ) ; 


GAL  BURGERQUICKDEQN 

function  da=gal  burgerquickdeqn (t, a, eps) 

%dif ferential  equation  referenced  in  GAL  BURGERDEQN 

global  P2  eps2 
NN=length (a) ; 
for  j  =1 : NN 

nonlin ( j ) =a ' *P2 ( : , : , j ) *a; 

end 

da=-pi A2*eps2* ( [1 :NN] ' . A2) . * a- 2* nonlin ' ; 


GAL  BURGERDISPLAY 

function  M=gal  burgerdisplay (t, x, u) 
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%utility  for  plotting  the  output  of  burger's  equation  files. 

O, 

o 

%inputs:  t--  times  at  which  solution  is  simulated,  equal  to  tvec  or 

else  output  of  ode45 

%  x —  points  along  x  axis  where  x  is  evaluated,  equal  to  xvec 

%  u —  output,  values  of  u  at  points  along  x  and  t.  u  has 

length (t)  rows  and  length (x)  columns 

O, 

o 

%Output :  M-  a  movie  that  shows  the  evolution  of  burger's  equation  over 

time 

%  figure (1)  contour  plots  of  u(x,t) 

%  figure (2)  surface  plot  of  u(x,t),  to  show  changes  over  time 

O, 

o 

%see  also:  GAL_BURGER,  GAL_BURGERDEQN ,  GAL_BURGERQU I CK , 

GAL_BURGERQU I CKDEQN 

cla 

%routine  for  the  titles 
blab=input ( ' what  basis?  \n','s'); 
tspan=input ( ' what  time  span?  \n','s'); 
xspan=input ( ' what  domain?  \n '  ,  's'); 
n=input('how  many  basis  functions?  \n') 

%contour  plots  of  u(x,t)  v.s.  x 

figure ( 1 ) 

view  (2) 

plot (x, u ' , ' k ' ) 

xlabel ( ' x ' ) 

ylabel ( ' u (x, t) ' ) 

title (' Burger '' s  equation  contour  plot') 

legend (['n=',int2str(n)], [ ' basis= ' , blab] , [' time= ' , tspan] , [ ' x= ' , xspan] ) 

%surface  plot  of  u(x,t)in  space  and  time 

figure (2 ) 

surf (x, t ' , u) 

xlabel ( ' X ' ) 

ylabel ( ' time ' ) 

z label ( ' u (x, t) ' ) 

title (' Burger '' s  equation  surface  plot') 

legend (['n=',int2str(n)], [ ' basis= ' , blab] , [' time= ' , tspan] , [ ' x= ' , xspan] ) 

j=input ( 'movie?  1  if  no,  2  if  yes') 
if  j==l 

return 

end 

%movie  routine 

figure ( 3 ) 

hold  off 

plot (x, u  (1,  : ) ) 

v=axis 

xlabel ( ' x ' ) 

ylabel ( ' u (x, t)  '  ) 

title (' Burger '' s  equation  movie') 
m=length (t) ; 

M=moviein (m) ; 
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for  j=l :m 

plot (x, u  ( j ,  : ) ) 
axis (v) 

M ( : ,  j ) =getframe; 
hold  off 

end 


GAL  HEAT 

%This  script  file  finds  a  solution  for  the  Heat  Equation  in  one 
dimension  using  the  Galerkin  Method 
%on  the  interval  [0,1] 

%and  then  plots  the  solution  at  various  times. 

%the  amount  of  time  is  an  input  t 

%the  number  and  type  of  basis  functions  are  inputs  n,  bn. 

%the  boundary  conditions  u(t,0)  and  u(t,l)  are  inputs  ut0,utl. 

clear  all 

%n  is  the  number  of  basis  functions 
n=4 

%bn  indicates  the  type  of  basis  function 
%1  gives  the  basis  cos (n*pi*x) 

%2  gives  the  basis  sin (n*pi*x) 

%3  gives  the  chebyshev  polynomial  basis  T(n) 
bn=2 

%these  are  the  boundary  conditions  u(t,0)=ut0  and  u(t,l)=utl 

ut0=l 

utl=2 

%%degree  of  PDE  for  t 
tdeg=l 

%t  is  the  time  we  let  the  simulation  run 
t=10 


%sets  up  a  as  a  matrix  of  the  functions  in  time  and  their  derivatives 
%also  sets  up  bv  as  a  vector  of  symbolic  basis  functions 
x=sym ( ' x ' ) ; 
for  j=0:n 

for  k=0:tdeg 

a ( k+1 , j  +1 ) =sym ( [ ' a ' , int2str ( j )  ,int2str(k)  ]  )  ; 

end 

bv ( j  +1 ) =basis ( j , x, bn) ; 
end 
a 

bv 

%gives  a  generic  linear  function  to  satisfy  the  boundary  conditions 
u(0,t),  u(l,t).  We  will  call  this 
%the  inital  condition  function  g=u(x,0) 
g= (utl-utO) *x+ut0 
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%plugs  the  assumed  solution  into  the  left  (eqnl)  and  right  (eqnr)  sides 
of  the  equation. 

%k  is  the  heat  constant 

%f  is  the  source  (can  be  in  the  symbolic  variable  x  if  desired) 


k=3 

f=l 

ya=a ( 1 , : ) . *bv; 
ya=  [g,  ya] 
dyadt=a ( 2 , : ) . *bv 

eqnr=sum ( k*dif f (ya, x, 2 ) ) +f 
eqnl=sum (dyadt) 


%multiplies  both  the  left  and  right  hand  sides  by  the  basis  function 
and  then  integrates  from  0  to  1 
for  j=0:n 

vl  (j+1) =int (eqnl . *bv (j+1)  ,  x,  0, 1)  ; 
vr(j+l)=int (eqnr*bv (j+1)  ,x,0,  1)  ; 

end 

global  vl 
global  vr 
vl 
vr 

%since 

bc=zeros  (n, 1 ) 

[t , avec] =ode45 ( ' gal_heatDEQN ' , 10 , be) ; 
xv=0 : .01:1; 

%avec ( 1 , : ) . *bv 
%sum(avec(l, :) .*bv) 

%subs (sum (avec ( j , : ) . *bv) , x, xv) 
cla 

hold  on 

for  j=l : length (t) 

plot (xv, xv+l+subs ( sum ( [ 0 , avec ( j , : ) ] . *bv) , x, xv) ) 

end 

%computes  value  for  differential  equation  using  the  inital  conditions 


GAL  HEATDEQN 

function  s=gal  heatDEQN (t, v) 

s (1) =2* (-1/ (2*pi) *3*piA3*v (1) +4/ (2*pi) ) ; 

s (2) =2*-6*piA2*v (2) ; 

s (3) =2*-l/ (6*pi) * (81*piA2*v(3) -4) ; 

s (4) =2*-24*piA2*v (4) ; 

s— s  '  ; 

%global  vl 
%global  vr 
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%n=length (vl) ; 

%dpolynorm=dpoly/dpoly ( 1 )  ; 

%f or  j  =1 : n 

%  scoef f=sym2poly (vl ( j ) )  ; 

%  vpoly=sym2poly (vr ( j ) )  ; 

%  while  length (vpoly<2 ) ; 

%  vpoly= [ 0 , vpoly]  ; 

%  end 

%  vcoef f=vpoly ( 1 )  ; 

%  vscalar=vpoly (2 )  ; 

%  s  ( j ) =vcoeff / scoeff *v ( j ) ivscalar/ scoeff; 

%end 

%s 


GAL  PEARSON2D 

function  [u, v, ab, X, t, init] =gal  pearson2D (NN, F2 , k2 , Xvec, tvec) 

% [u, v, ab, X, t, init] =gal  pearson2D (NN, F2 , k2 , Xvec, tvec) 

%[u,v,ab,X,t,init] =gal_pearson2D (4,  .0152,  .04,  [  [-1 : . 01 : 1 ]  ',  [  - 

1:  .01:1]  ' ] ,  1000) ; 

%this  function  solves  Pearson's  2  morphogen  reaction-diffusion  equation 
in  2  dimensions 

%using  the  galerkin  spectral  method  with  periodic  boundary  conditions. 
The  resulting  ODE's  are  solved 

%using  ODE45  (runge  kutta  4  method,  forward  time  differencing) . 

O, 

o 

%the  equations  are  defined  as 

o, 

o 

%ut=du*laplacian (u) -uvA2+F (1-u) 

%vt=dv*laplacian (v) +uvA2- (F+k)  v 

O, 

o 

%u,v:  The  concentrations  of  each  morphogen 

%F:  Dimensionless  feed  rate 

%k:  Dimensionless  rate  constant 

%du,dv:  Diffusion  constants  for  u  and  v  respectively 

O, 

o 

%%Ref:  Complex  Patterns  in  a  Simple  System,  John  E.  Pearson,  Science, 

V261,  5118,  9  Jul  1993  189-192 

o, 

o 

%Inputs : 

%NN2=number  of  basis  functions  used 
%F2 , k2-scalar  parameters 

%Xvec-  2  column  space  vector,  column  one  is  the  X  coordinate,  column  2 
is  the  Y  coordinate 

%tvec-  time  vector  discretizing  time  for  ODE45 

O, 

o 

%Outputs : 

%u,v-  Morphogen  concentrations 
%X-  Xvec 
%t-  tvec 

o, 

o 

%See  GAL_PEARSON2DDEQN 

%initialize  global  variables  F, k, du, dv, parameters,  NN,  basis  number, 

M2,  sparse 
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%inner  product  matrix 

global  F  kk  du  dv 

F=F2 ; 

kk=k2 ; 

du= . 00002 ; 

dv=. 00001; 


%finds  the  vector  of  constants  for  the  Laplacian,  G 

global  G 

G=-NN:NN; 

G=G. *2; 
n=l; 

for  k=G; 

G2 (:,n)=G'+k; 
n=n+l ; 

end 

G=G2  (  :  )  ; 

%finds  initial  condition 
tic, 

for  p=-NN : NN 

for  q=0 : NN  %-NN : NN 

%initu (p+NN+1, q+NN+1) =quadl (@initcondu, - 
1,1,  [],  [] ,p) *quadl (@initcondu, -1,1,  []  ,  [ ]  ,  q)  ; 

initu (p+NN+1 , q+1 ) =dblquad ( @initcondu, -1,1, -1,1, [], @quadl ,  p,  q- 

NN)  ; 

initu (NN+l-p, 2*NN+l-q) =conj (initu (p+NN+1, q+1) ) ; 

%initv (p+NN+1, q+NN+1) =quadl (@initcondv, - 
1,1,  [],  [] ,p) *quadl (@initcondv,  -1,1,  []  ,  []  ,  q)  ; 

ini tv (p+NN+1 , q+1 ) =dblquad ( @initcondv, -1,1, -1,1, [], @quadl , p, q- 

NN)  ; 

initv (NN+l-p, 2*NN+l-q) =conj  (initv (p+NN+1,  q+1)  )  ; 

end 

end, 

initu=initu ( : ) / 4 ; 
initv=initv ( : ) / 4 ; 
init= [initu; initv] ; 

' init  matrix  finished', 
toe 

%Calling  the  ODE  solver 

% [t, ab] =ode45 ( ' gal_pearson2Dexpdeqn ' , tvec,  init) ; , 
tic, 

[t, ab] =ode45 ( ' gal_pearson2Ddeqn '  ,  tvec,  init)  ; 

' ode '  ' s  solved '  , 
toe 

%constructing  the  solution 
tic, 

x=Xvec  ( : , 1 )  '  ; 
y=Xvec  (  : , 2 )  ; 

X=Xvec; 

u=zeros (length (x) , length (y) , length (t) ) ; 
v=zeros (length (x) , length (y) , length (t) ) ; 
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for  j=l : length (t) 

for  k=l : size (ab, 2 ) /2 
p=mod (k, 2*NN+1) ; 
if  p==0 

p=2*NN+l; 

end 

p=p- (NN+1 ) ; 
q=ceil (k/ (2*NN+1) ) ; 
q=q- (NN+1 ) ; 
modey=exp (i*pi*q*y)  ; 
modex=exp ( i*pi*p*x) ; 

u ( : , : , j ) =u ( : , : , j ) +ab ( j , k) *modey*modex; 
v ( : ,  : , j ) =v ( : ,  : , j ) +ab ( j , k+size (ab, 2) / 2)  *modey*modex; 

end 

end, 

'solution  constructed', 
toe 


function  z=initcondu (x, y, p, q)  %,k) 

%z=exp (~i*pi*k*x)  .* (l+.l* (x.A2-l)  )  ; 

%z=l- . 5* (exp (-20* (x-.l)  ,A2) ) * (exp (-20* (y) A2) ) ;  %-l/16*exp (-20* (x- 

. 1) . A2) *exp (-30* (y-.l) A2) ; 

z=l- . 5* (exp (-20* (x+.7) . A2) ) * (exp (-20* (y-.75) A2) ) ;  %-l/16*exp (-20* (x- 

.  1)  .  A2) *exp (-30* (y-.l) A2)  ; 
z=exp (-i*pi* (p*x+q*y) ) .  *z; 

function  z=initcondv (x, y, p, q)  %,k) 

%z=exp (-i*pi*k*x)  .  * (.1* (l-x.A2)  )  ; 

%z= . 25* (exp (-20* (x- . 05 ) . A2 ) ) * (exp (-20*yA2) ) ;  %+l/16*exp (-20* (x- 

.15) . A2) *exp (-30* (y-.l) A2) ; 

z= . 25* (exp (-20* (x+.7) ,A2) ) * (exp (-20* (y-.65) A2) ) ;  %+l/16*exp (-20* (x- 

.15) . A2) *exp (-30* (y-.l) A2) ; 
z=exp (-i*pi* (p. *x+q*y) ) . *z; 


GAL  PEARSON2DDEQN 

function  dabdt=gal  pearson2Ddeqn (t, ab) 

%two  dimensional  head  equation  differential  equation  called  in 
gal  heat2D. 

%this  program  sets  up  the  system  of  differential  equations  inolving  the 
a  (m,  n)  (t) 

global  F  kk  G  du  dv 

NN= (sqrt (length (ab) / 2) -1) / 2; 

a=ab ( 1 : length (ab) / 2)  ; 
b=ab (length (ab) /2  +  1 :end)  ; 

%does  the  a's  without  the  nonlinear  part.  The  extra  4*Fv  is  added 
because  it  is  the  (phi(0,0),F)  term.  For  every 
%other  equation,  this  disappears 
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Fv=zeros (size (a) ) ; 

Fv ( ( (2*NN+1 ) A2+l ) /2 ) =F; 


dabdtu= (-du*piA2*G-F) ,*a+Fv;  %-sum (nonlinear ) ; 

%does  the  b's 

dabdtv= (-dv*piA2*G-F-kk) ,*b;  %+sum (nonlinear ) ; 
for  j  =1 : 2  *NN+1 

bm( : , j ) =b ( ( j-1) * (2*NN+1 ) +1 : j * (2*NN+1) )  ; 

end 

bm2= [zeros (2*NN, 6*NN+1 ) ;  [zeros (2  *NN+1 , 2  *NN) , bm, zeros (2*NN+1 , 2  *NN) ] ; zero 
s  (2*NN, 6*NN+1 ) ] ; 

for  q=-NN : NN 

for  p=-NN : NN 

for  k2=-NN : NN 

for  kl=-NN : NN 

for  m2=-NN : NN 

for  ml=-NN : NN 

bmpq ( kl+NN+l+ ( k2+NN) * (2*NN+1) ,ml+NN+l+ (m2+NN) * (2*NN+1) ) =bm2 ( (p-kl- 
ml)+3*NN+l, (q-k2-m2) +3*NN+1) ; 

%bmpq ( kl k2 , mlm2 ) 

*bm2 (p-kl-ml , q-k2-m2 ) 

end 

end 

end 

end 

nonlin (p+NN+l+ (q+NN) * (2*NN+1 ) ) =sum ( sum ( (a*b . ' ) . *bmpq) ) ; 
%nonlin(p,q)  =  sum(akl,k2  x  bml*m2  x  b(p-ml- 

kl ) , (q-m2-k2 ) 
end 

end 

dabdtu=dabdtu-nonlin  ( : ) ; 
dabdtv=dabdtv+nonlin  (:); 
dabdt= [dabdtu; dabdtv] ; 

GAL  1DPDEDISPLAY 

function  M=gal  lDpdedisplay (t, x, u) 

%type  gal_lDpdedisplay (t, x, u)  to  dislpay  the  results  as  a 

%contour  map,  a  surface  plot  or  a  movie 

%plotting 

cla 

%routine  for  the  titles 
tit=input ( ' input  Title\n ' ,  ’s'); 
blab=input ( ' what  basis?  \n','s'); 
tspan=input ( ' what  time  span?  \n’,’s'); 
xspan=input ( ' what  domain?  \n '  ,  ’s’); 
n=input('how  many  basis  functions') 

%contour  plots  of  u(x,t)  v.s.  x 
figure ( 1 ) 
view  (2) 
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plot (x, u ' ) 
xlabel ( ' x ' ) 
ylabel ( ' u (x, t) ' ) 
title (tit) 

legend ( [ ' n= ' , int2str (n) ] , [ ' basis= ' , blab] , [' time= ' , tspan] , [ ' x= ' , xspan] ) 

%surface  plot  of  u(x,t)in  space  and  time 

figure (2 ) 

surf (x, t ' , u) 

xlabel ( ' X ' ) 

ylabel ( ' time ' ) 

z label ( ' u (x, t) ' ) 

title ( [tit, ' surface  plot']) 

legend (['n=',int2str(n)], [ ' basis= ' , blab] , [' time= ' , tspan] , [ ' x= ' , xspan] ) 

j=input ( 'movie?  1  if  no,  2  if  yes') 
if  j==l 

return 

end 

%movie  routine 

figure ( 3 ) 

hold  off 

plot (x, u  (1,  : ) ) 

v=axis 

xlabel ( ' x ' ) 

ylabel ( ' u (x, t)  '  ) 

title ( [tit, '  movie ' ] ) 

m=length  (t) ; 

M=moviein (m) ; 

for  j=l :m 

plot (x, u ( j , : ) ) 
axis (v) 

M ( : , j ) =getframe; 
hold  off 

end 


GAL  2DPDEDISPLAY 

function  M=gal  2Dpdedisplay (u, X, t) 

%produces  movie  images  for  given  2D  pde  data  in  the  form: 
%t=  time  vector 

%X—  space,  x  in  column  1,  y  in  column  2 

%u=  u(x,y,t)=  the  value  of  the  function  at  x,y  and  t 

cla 

tit=input ( ' Input  Title\n ' , 's') 

%movie  routine 

f=input (' Surface (2 )  or  Contour  ( 1 ) \n ' ) 
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if  f==l 

hold  off 
x=X ( : , 1 )  ; 
y=X ( : , 2 )  ; 
xlabel ( ' x ' ) 
ylabel ( ' y ' ) 
title (tit) 
m=length  (t) ; 

M=moviein (m) ; 
for  j  =1 : m 

contour (x, y, u ( :  ,  :  ,  j  )  ) 
title (tit) 

M ( : , j ) =getframe; 
hold  off 

end 


elseif  f==2 

hold  off 
x=X ( : , 1 )  ; 
y=X ( : , 2 )  ; 
u=real (u) ; 

v=  [  -1 ,  1 ,  -1 ,  1 ,  min  (min  (min  (u)  )  )  ,  max  (max  (max  (u)  )  )  ]  ; 

xlabel ( ' x ' ) 

ylabel ( ' y ' ) 

z label ( ' u (x, y, t) ') 

title (tit) 

m=length (t) ; 

M=moviein (m) ; 
for  j=l :m 

mesh (x, y, u ( : , : , j ) ) 
axis (v) 
title (tit) 

M ( ; , j ) =getframe; 
hold  off 

end 

else 

'wrong  input' 


end 


GAL  ABRECONSTRUCT 

function  [u, v] =gal_abreconstruct (ab, t, init) 

%this  function  reconstructs  u  and  v  from  given  ab  vector,  t,  and  init 

values  from 

%2D  pde  results: 

o, 

o 

%inputs:  ab-  length (t) xlength (init)  vector  of  coefficient  functions 

%  t-  times  at  which  coefficient  functions  are  evaluated 

%  init-  initial  conditions  for  the  ab  same  as  ab ( : , 1 ) 
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% syntax :  [u, v] =gal_abreconstruct (ab, t, init) 

tic, 

x= [ - 1 : .01:1]  ; 
y=x '  ; 


u=zeros (length (x) , length (y) , length (t) ) ; 
v=zeros (length (x) , length (y) , length (t) )  ; 
NN= (sqrt (length (ab) / 2) — 1) /2; 


for  j=l : length (t) 

for  k=l : size (ab, 2 ) /2 
p=mod (k, 2*NN+1) ; 
if  p==0 

p=2*NN+l; 

end 

p=p- (NN+1 ) ; 
q=ceil (k/ (2*NN+1) ) ; 
q=q- (NN+1 ) ; 

modey=exp (i*pi*q*y)  ; 
modex=exp (i*pi*p*x) ; 

u ( : ,  : , j ) =u  ( : ,  : , j ) +ab ( j , k) *modey*modex; 
v ( : ,  : , j ) =v ( : ,  : ,  j ) tab ( j , k+size (ab, 2) / 2)  *modey*modex; 

end 

end, 

'solution  constructed', 
toe 


Gray-Scott  Analysis 

GAL2DGSN0EIGENVALS 

%eigenvalue  movie  for  the  diffusionless  Grey  Scott  Model 

o, 

o 

%  a_t=-a*bA2-F (1-a) 

%  b  t=a*bA2+ (F+k) b 

o, 

o 

%This  routine  creates  a  movie  of  the  eigenvalues  of  the  derivative 
matrix  of  { a_t (a, b) , b_t (a, b) } 

%evaluated  at  the  equilibrium  points  of  the  system  in  the  complex  plane 
when  k  is  held  constant 

%at  .04  and  F  varies  in  order  to  examine  local  stability  of  the 
equilibrium  points. 

cla 

hold  on 
k= .  04 

xlabel ( ' real  axis ' ) 
ylabel (' imaginary  axis') 
n=l 

Fi=. 01 


Ff=. 17 
df=. 0001 
for  j=Fi:df:Ff 
F=j; 

eq0= [ -F, 0 ; 0 ,  - (F+k)  ]  ; 
eig0=eig (eqO) ; 
figure  ( 1 ) 

axis  (  [-.25,  .15, -.15,  .15] )  ; 
hold  on 

plot (eigO, zeros (size (eigO) ) , ' x ' ) 
r=roots ( [F+k, -F, FA2+F*k] ) ; 
r=sort ( r ) ; 
b=r  (1)  ; 
a= (F+k) /b; 

eql=[-bA2-F, -2*a*b;bA2, 2*a*b- (F+k) ]  ; 
eigl=eig (eql )  ; 
figure (2 ) 

axis (  [-.25,  .15, -.15,  .15] )  ; 
hold  on 

if  isreal (b) ==1 

if  isreal (eigl ) ==1 

plot (eigl, zeros (size (eigl) ) , 'x' ) 

else 

plot (eigl , ' x ' ) 

end 

end 

b=r  (2)  ; 
a= (F+k) /b; 

eq2=[-bA2-F, -2*a*b;bA2, 2*a*b- (F+k) ]  ; 
eig2=eig (eq2 ) ; 
figure ( 3 ) 

axis (  [-.25,  .15, -.15,  .15] )  ; 
hold  on 

if  isreal (b) ==1 

if  isreal (eig2 ) ==1 

plot (eig2, zeros (size (eig2) ) , 'x' ) 

else 

plot (eig2 , ' x ' ) 

end 

end 

%M ( : , n) =getf rame; 
n=n+l ; 

end 

figure  ( 1 ) 

plot  ([0,0],  [-.15,  .15],  ' k ' ) 
plot ( [- . 25, .15],  [0,0],  ' k ' ) 
xlabel ( ' real  axis ' ) 
ylabel (' imaginary  axis') 

text (-. 25, . 15, ' Eigenvalues  of  the  stability  matrix  of  (a  0,b  0)') 
text (-.23, .125, [ ' k= ' , num2str ( k ) , ' ,  F  increases  from  ' , num2str (Fi) , 
' , num2str (Ff ) ] ) ; 

figure (2 ) 

plot  ([0,0],  [-.15,  .15],  ' k ' ) 
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plot (  [-.25,  . 15] ,  [0, 0] ,  '  k'  ) 
xlabel ( ' real  axis ' ) 
ylabel (' imaginary  axis') 

text (-. 25, . 15, ' Eigenvalues  of  the  stability  matrix  of  (a  l,b  1)') 
text (-.23, .125, [ ' k= ' , num2str (k) , ' ,  F  increases  from  ' , num2str (Fi) , '  to 
' , num2str (Ff ) ] ) ; 

figure ( 3 ) 

plot  ([0,0],  [-.15,  .15],  ' k ' ) 
plot  (  [- . 25,  .15],  [0,0],  ' k ' ) 
xlabel ( ' real  axis ' ) 
ylabel (' imaginary  axis') 

text (-. 25, . 15, ' Eigenvalues  of  the  stability  matrix  of  (a_2,b  2)') 
text (-.23, .125, [ ' k= ' , num2str (k) , ' ,  F  increases  from  ' , num2str (Fi) , '  to 
' , num2str (Ff ) ] ) ; 


GAL2DGSHOPFFIND 

%find  value  of  F  where  eig2  crosses  the  axis  and  becomes  stable 
k=.04-->F=. 0152786 

%answer  gives  estimate  that  is  less  than  the  true  answer  by  no  more 

than  tol 

cla 

hold  on 
n=l 

tol=10 A-10 

for  k=0 : .0005: .0625 

%finds  starting  F  value 
F=roots ( [-4, 1-8 *k, -4*kA2, 0] ) ; 

F=F (F>0) ; 

Fmin=min (F) ; 

Fmax=max (F) ; 

eq3bot (n, : ) = [k, Fmin] ; 

eq3top  (n,  : )  =  [k, Fmax] ; 

%adds  a  bit  to  starting  value  to  be  sure  that  the  equilibrium  is 
real  but  not 

%enough  to  make  it  stable 
F=Fmin+. 00001; 

%initializes  eig2  for  while  loop 
r=roots ( [F+k, -F, FA2+F*k] ) ; 
r=sort ( r ) ; 
b=r  (2)  ; 
a= (F+k) /b; 

eq2=[-bA2-F, -2*a*b;bA2, 2*a*b- (F+k) ]  ; 
eig2=eig (eq2 )  ; 

%set  initial  step  value 
d= . 0  0 1 ; 

%outer  while  loop  decreases  step  value  until  it  is  within  tolerance 
while  d>tol 

%inner  while  loop  tests  eig2  for  stability 
while  real(eig2)>0 
F=F+d; 

r=roots ( [F+k, -F, FA2+F*k]  )  ; 
r=sort ( r ) ; 
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b=r (2) ; 
a= (F+k) /b; 

eq2=[-bA2-F,-2*a*b;bA2,2*a*b-(F+k)  ]  ; 
eig2=eig (eq2 )  ; 

end 

%sets  F  to  just  above  where  it  broke  the  while  loop  so  it  can 
go  through  again  at 

%higher  tolerance  with  a  reinitialized  eig(2) 

F=F-d; 

r=roots ( [F+k, -F, FA2+F*k]  )  ; 
r=sort (r) ; 

b=r  (2) ; 
a= (F+k) /b; 

eq2=[-bA2-F, -2*a*b;bA2, 2*a*b- (F+k) ] ; 
eig2=eig (eq2 ) ; 

%decreases  step  value 
d=d* . 1; 

end 

%stores  answers 
Fhopf (n) =F; 
n=n+l ; 
end 

eq3= [eq3bot ( [2 : end] ,  :  )  ; f lipud (eq3top) ] ; 
plot (eq3 {:,!), eq3 ( :  ,  2) ) 
kvec=eq3bot (2 : end, 1 ) ; 

Fhopf=Fhopf (2 : end) ; 
plot ( kvec ' , Fhopf, ' -- ' ) 


GAL2DGSFIELDPLOTS 

%plots  fieldlines  and  equilibrium  points.  Can  be  altered  to  make  a 

movie  as  F  is  varied 

n=l 

figure (2 ) 
for  j=.0154 
cla 

hold  on 

F=j 

k= .  04 

umin=0 

umax=l 

vmin=0 

vmax=l 

spa=20 

time=100 

%plotting  the  roots  of  the  equation 
plot  (1, 0,  ' ro ' ,  ' markers ize ' , 15) 
eq0=[-F,0;0,-(F+k) ] 
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eigO=eig (eqO ) 

r=roots ( [F+k, -F, FA2+F*k] ) ; 
r=sort ( r ) ; 
if  isreal(r)==l 

plot( (F+k) . /r ( 1 ) ,r(l) , ' r+ ' , ' markers ize ' , 15 ) 
plot ( (F+k) . /r (2 ) , r (2 ) , ' ro ' , ' markers ize ' , 15 ) 

end 

text (umax- . 1* (umax-umin) , vmax- . 1* (vmax- 
vmin) , [num2str (umax) , ' , ' , num2str (vmax) ] ) ; 
text (umin+ . 1* (umax-umin) , vmin+ . 1* (vmax- 
vmin) , [num2str (umin) , ' , ' , num2str (vmin) ] ) ; 


b=r  (1)  ; 
a- (F+k) /b; 

eql=[-bA2-F, -2*a*b;bA2, 2*a*b- (F+k) ] 
eigl=eig (eql ) 

b=r  (2)  ; 
a= (F+k) /b; 

eq2=[-bA2-F, -2*a*b;bA2, 2*a*b- (F+k) ] 
eig2=eig (eq2 ) 

for  j =linspace (umin, umax, spa) 

for  m=linspace (vmin, vmax, spa) 

[t, u] =ode45 ( ' gal_pearsonOdeqn ' , time, [ j , m] , [ ] , F, k) ; 
plot (u ( : , 1) , u ( : , 2) ) 

end 

end 

xlabel ( ' u ' ) 
ylabel ( ' v ' ) 

title ( [ ' F= ' , num2str (F) , '  k= ' , num2str (k) ,  '  time= ' , num2str (time) ] ) 

axis ( [umin, umax, vmin, vmax] ) 

%M ( : , n) =getframe; 

n=n+l ; 

end 


GALPE  ARSON  ODE  QN 

%function  used  in  Gray  Scott  Analysis  functions  for  the  N=0 
differential  equation 

function  dabdt=gal_pearsonOdeqn (t, ab, flag, F, k) 
dabdt (1,1) =-ab ( 1 ) . *ab (2) . A2+F* (1-ab (1) ) ; 
dabdt (2,1) =ab ( 1 ) . *ab (2 ) . A2- (F+k) *ab (2 ) ; 
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GAL  DDINSTABILITY 

clear  all 
F= . 0154 
k= .  04 
du=. 00002 
dv=. 00001 

n=l  ; 

modes=0 : 1 0 : 5000 ; 
for  mn=modes; 

b= (F+sqrt (FA2-4*F* (F+k) A2) ) /  (2* (F+k)  )  ; 
a= (F+k) /b; 

M= [-bA2-F-mn*du, -2*a*b;bA2, 2*a*b- (F+k) -mn*dv] ; 
eva=eig (M) ; 
y (n) =real (eva (2 ) ) ; 
n=n+l ; 

end 

plot (modes , y) ; 


